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Abstract. We prove convergence of discrete duality finite volume (DDFV) 
schemes on distorted meshes for a class of simplified macroscopic bidomain 
models of the electrical activity in the heart. Both time-implicit and linearised 
time-implicit schemes are treated. A short description is given of the 3D DDFV 
meshes and of some of the associated discrete calculus tools. Several numerical 
tests are presented. 
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1. Introduction 



We consider the heart of a living organism that occupies a fixed domain f2, 
which is assumed to be a bounded open subset of K 3 with Lipschitz boundary 
9f2. A prototype model for the cardiac electrical activity is the following nonlinear 
reaction-diffusion system 



where Q denotes the time-space cylinder (0, T) x f2. 

This model, called the bidomain model, was first proposed in the late 1970s by 
Tung [M] and is now the generally accepted model of electrical behaviour of cardiac 
tissue (see Henriquez [33], Keener and Sneyd [ID]). The functions Ui — Ui(t,x) 
and u e = u e (t,x) represent the intracellular and extracellular electrical potentials, 
respectively, at time t £ (0,T) and location x £ £1. The difference v = m — u e 
is known as the transmembrane potential. The conductivity properties of the two 
media are modelled by anisotropic, heterogeneous tensors Mj(x) and M e (x). The 
surface capacitance of the membrane is usually represented by a positive constant 
c. m ; upon rescaling, we can assume c m = 1. The stimulation currents applied to 
the intra- and extracellular spaces are represented by an L 2 (Q) function / app = 
Ia,pp{t,x). Finally, the transmembrane ionic current h[v] is computed from the 
potential v. The system is closed by choosing a relation that links h[v] to v and 
specifying appropriate initial-boundary conditions. We stress that realistic models 
include a system of ODEs for computing the ionic current as a function of the 
transmembrane potential and a series of additional "gating variables" aiming to 
model the ionic transfer across the cell membrane (see, e.g., S3 S3 HO]) ■ This 
makes the relation h = h[v] non-local in time. 

Herein we focus on the issue of discretisation in space of the bidomain model. The 
presence of the ODEs, some of them being quite stiff, greatly complicates the issue of 
discretisation in time. It also results in a huge gap between theoretical convergence 
results and the practical computation of a reliable solution. We surmise that the 
precise form of the relations that link h[v] to v is not essential for the validation of 
the space discretisation techniques. Therefore, as in [5] [3] , we study ([I]) under the 
greatly simplifying assumption that the ionic current is represented locally, in time 
and space, by a nonlinear function h(v). However, such a simplification allows to 
mimic, to a certain extent, the depolarisation sequence in the cardiac tissue, taking 
the ionic current term h[v] to be a cubic polynomial (bistable equation); this choice 
models the fast inward sodium current that initiates depolarisation (cf., e.g., |17j). 

In the context of electro-cardiology the relevant boundary condition would be 
a Neumann condition for the fluxes associated with the intra- and extracellular 
electrical potentials: 



It serves to couple the heart electrical activity with the much weaker electrical 
phenomena taking place in the torso. The simplest case is the one of the isolated 
heart, namely S; !e = 0. For the mathematical study we are heading to, we consider 
rather general mixed Dirichlet-Neumann boundary conditions of the form 

(2) u l:C — g i>e on (0,T) x r fl , M, iE (j;)Vvn = V on(0 1 T)xT N , 

where dfl is partitioned into sufficiently regular parts and Tp, and n denotes 
the H 2 -a.e. defined exterior unit normal vector to the Neumann part Tjy of the 
boundary d£l. To keep the analysis simple, let us assume that Sj, e £ L 2 ((0, T) xfjv); 



(1) 




{t,x) £ Q, 
(t,x) £ Q, 



M iye (x)\/ui, e ■ n = Sj,, 



i,e 



on (0,T) x dfl. 
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for g^ e , we assume g^ e € L 2 (0,T; H 1 / 2 (T £>)) (in fact, we consider g^ e extended to 
^(O/T;^ 1 ^)) functions). 

Regarding the initial data, we prescribe only the transmembrane potential: 

(3) v(Q,x) — vq{x), x £ £1. 

Clearly, ([lj and ^ are invariant under the simultaneous change of Ui^ e into 
Ui ie + k, k e M. In the case dQ = Tjv, = 0, also ([2| is invariant under this 
change; therefore, for the sake of being definite, we normalise u e by assuming 

(4) whenever T D = 0, / u e (t, ■) = for a.e. t € (0,T). 

Jo 

It is easy to see that the existence of solutions to |l]),([3| requires the compatibility 
condition 

(5) whenever T D = 0, / Si(t, •) + s e (t, •) = for a.e. i € (0, T). 

Jon 

Notice that the diffusion operators Mi )e (x)Vui ie in ([I]) are linear in the gradient 
Vu^ei heterogeneous & anisotropic, and time-independent; these assumptions seem 
to be sufficiently general to capture the phenomena of the electrical activity in the 
heart. More general models with time-dependent and nonlinear in Vui je diffusion 
of the Leray-Lions type were studied in [8]. Here we assume that (M. iiB (x)} Q is 
a family of symmetric matrices, uniformly bounded and positive definite: 

3 7 for a.e. x e ft, V£ G K 3 , -|£| 2 < (M i , e (a;)0 • £ < 7 |£| 2 . 

7 

In particular, we have Mj )6 € L°°(f2). 

Now let us describe in detail the ionic current function h = h(v). We assume 
that h : M — > K is a continuous function, and that there exist r G (2, +oo) and 
constants a,L,l > such that 

(6) -\v\ r <\h(v)v\<a(\v\ r + 1), 

a 

(7) h : z i — y h{z) + Lz + I is strictly increasing on K, with lim h(z)/z = 0. 
For the later use, we set 

b : z n- 6(0) = 0. 

It is rather natural, although not necessary, to require in addition that 
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(8) Vz.seR (h(z) - h(s))(z - s) > ^(l + \z\ + \s\) r - 2 \z-s\ 2 . 

According to [THJ [TS] , the most appropriate value is r = 4, which means that the 
non-linearity h is of cubic growth at infinity. Assumptions (|6]),([7| are automatically 
satisfied by any cubic polynomial h with positive leading coefficient. 

A number of works have been devoted to the theoretical and numerical study of 
the above bidomain model. Colli Franzone and Savare [19] prove the existence of 
weak solutions for the model with an ionic current term driven by a single ODE, by 
applying the theory of evolution variational inequalities in Hilbert spaces. Sanfelici 
[5T] considered the same approach to prove the convergence of Galerkin approxi- 
mations for the bidomain model. Veneroni in |55j extended this technique to prove 
existence and uniqueness results for more sophisticated ionic models. Bourgault, 
Coudiere and Pierre [T3] prove existence and uniqueness results for the bidomain 
equations, including the FitzHugh-Nagumo and Aliev-Panfilov models, by applying 
a semigroup approach and also by using the Faedo-Galerkin method and compact- 
ness techniques. Recently, Bendahmane and Karlsen [S] proved the existence and 
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uniqueness for a nonlinear version of the simplified bidomain equations ([!]) by using 
a uniformly parabolic regularisation of the system and the Faedo-Galerkin method. 

Regarding finite volume (FV) schemes for cardiac problems, a first approach is 
given in Harrild and Henriquez |31| , Coudiere and Pierre [22] prove convergence 
of an implicit FV approximation to the monodomain equations. We mention also 
the work of Coudiere, Pierre and Turpault [23j on the well-posedness and testing of 
the DDFV method for the bidomain model. Bendahmane and Karlsen [9] analyse 
a FV method for the bidomain model with Dirichlet boundary conditions, supply- 
ing various existence, uniqueness and convergence results. Finally, Bendahmane, 
Burger and Ruiz [10] analyse a parabolic-elliptic system with Neumann boundary 
conditions, adapting the approach in [S]; they also provide numerical experiments. 

In this paper, as in [9], we use a finite volume approach for the space discretisation 
of (JlJ and the backward Euler scheme in time. Due to a different choice of the finite 
volume discretisation, we drop the restrictions on the mesh and on the isotropic and 
homogeneous structure of the tensors M^ e imposed in [5]. We also consider general 
boundary conditions ([2| . The space discretisation strategy we use is essentially the 
one described and implemented by Pierre [35] and Coudiere et al. [231 [21] ■ More 
precisely, we utilise different types of DDFV discretisations of the 3D diffusion 
operator; in addition to the scheme of [49[[24], we examine the schemes described 
in [3 [351 12] (see also [I]) and [3T]. It should be noticed that 2D bidomain simulations 
on slices of the 3D heart are also of interest. The standard 2D DDFV construction 
can be applied to problem ([!]), ^ , ^ on 2D polygonal domains; the 3D convergence 
results readily extend to the 2D case. 

The DDFV approximations were designed specifically for anisotropic and/or non- 
linear diffusion problems, and they work on rather general (eventually, distorted, 
non-conformal and locally refined) meshes. We refer to Hermeline [3JJ ES3 [SHI [371 
135] . Domelevo and Omnes [55], Delcourte, Domelevo and Omnes [2S], Andreianov, 
Boyer and Hubert [5J, and Herbin and Hubert |32| for background information on 
DDFV methods. Most of these works treat 2D linear anisotropic, heterogeneous 
diffusion problems, while the case of discontinuous diffusion operators have been 
treated by Boyer and Hubert in [TS] . Hermeline [37J [35] treats the analogous 3D 
problems, [55] treats the Stokes problem, and the work [5] is devoted to the non- 
linear Leray-Lions framework. 

A number of numerical simulations of the full bidomain system (the PDE (nj) for 
Ui te plus ODEs for h[v]) coupled with the torso can be found in [43 ] 144" ! |23" 1 152 | [53] . 

Our study can be considered as a theoretical and numerical validation of the 
DDFV discretisation strategy for the bidomain model. For both a fully time-implicit 
scheme and a linearised time-implicit scheme, we prove convergence of different 
DDFV discretisations to the unique solution of the bidomain model ([I]). Then 
numerical experiments are reported to document some of the features of the DDFV 
space discretisations. A rescaled version of model ([I]), together with a cubic shape 
for v i — y h[v], is used to simulate the propagation of excitation potential waves 
in an anisotropic medium. In our tests, we combine 2D and 3D DDFV schemes 
for the diffusion terms with fully explicit discretisation of the ionic current term; 
thus numerical experiments validate this scheme, although we were not able to 
justify its convergence theoretically. Convergence of the numerical solutions towards 
the continuous one is measured in three different ways: the first two ones are 
aimed at physiological applications (convergence for the activation time and for 
the propagation velocity), whereas the third one corresponds to the norm used in 
Theorem |4.5| Implementation is detailed. Due to a large number of unknowns and a 
relatively large stencil of the 3D DDFV schemes, a careful preconditioning is needed 
for the bidomain system matrix that has to be inverted at each time step. The 
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preconditioning strategy we adopted here is developed in |50j : it provides an almost 
linear complexity with respect to the matrix size for the system matrix inversion. 
The preconditioning combines the idea of hierarchical matrices decomposition [111 
112] with heuristics referred to as the monodomain approximation |18j . 

The remaining part of this paper is organised as follows: In Section[2]we give the 
definition of a weak solution to ([I]) , ([2]) , (J3| . Moreover, we recast the problem into 
a variational form, from which we deduce an existence and uniqueness result. In 
Section[3]we describe one of the 3D DDFV schemes, while in Scction[2]we formulate 
two "backward Euler in time" & "DDFV in space" finite volume schemes, and state 
the main convergence results. The proofs of these results are postponed to Section 
[6j their basis being Section [5j where we recall some mathematical tools for studying 
DDFV schemes. Finally, Section [7] is devoted to numerical examples. 

2. Solution framework and well-posedness 

We introduce the space 

V = closure of the set {v € C°°(IR 3 ), v\ Fd = 0} in the iJ^fi) norm. 

In the case Yd = 0, we also use the quotient space Vq := V/{v G V, v = Const}. 
The dual of V is denoted by V', with a corresponding duality pairing (•, •). 
We assume that the Dirichlet data gi. e in |2]) are sufficiently regular, so that 

<7j e are the traces on (0, T) x of a couple of L 2 (0, T; H 1 ^)) functions 

(we keep the same notation for the functions gi^ e and their traces). For the sake of 
simplicity, we assume that 

the Neumann data s^ e belong to L 2 ((0,T) x Fn). 

Finally, we require that 

the initial function vq belongs to L 2 (Q). 

Definition 2.1. A weak solution to Problem ([!]), ([2| , ([3| is a triple of functions 

(9) (u h u e ,v) -M^R 3 s.t. u iie -g iie G L 2 (0,T;V), v = u t -u e , v G L r (Q), 

and such that |l]),([2|,([3| are satisfied in 2?'([0,T) x (fiur^v))- In the case = 0, 
we normalise u e by requiring Q. 



Remark 2.2. It is not difficult to show that Definition 2.1 is equivalent to a 
"variational" formulation of Problem p" ) , p| , ([3]) , in the spirit of Alt and Luckhaus 
PQ. Indeed, a triple (ui, u e , v) satisfying (9| is a weak solution of Problem ([!]), ^ , ^ 
if and only if , ^ , ^ are satisfied in the space L 2 (0, T; V') + L r '(Q). This means 
precisely that the distributional derivative dtv can be identified with an element of 
L 2 (0, T; V) + L r (Q), and with this identification there holds 



(10) 



f (d t v,tp}+ (( (M.^jVu, • V<p + h(v)<p) - I I Si <p= ft I app <p, 

JJQ Jo JT N JJQ 

/ (d t v,ip)- // (M e (ar) Vu e • Vip + h(v)ip) - / / s e tp= I^ip 

Jo JJQ Jo Jr N JJQ 



iQ 

for all ip e L 2 (0, T; V ) D U (Q) , and 

(dtv,<p)=- // vd t (p- / v (')<p(0,-) 
JJq Jn 

for all ip G L 2 (0, T, V) such that d t <p G L°°{Q) and (p(T, ■) = 0. 
We have the following chain rule: 
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Lemma 2.3. Assume thatv G L 2 (0, T; V)DL r (Q) and d t v G L 2 (0,T; V')+L r '(Q). 
Then 

T (d t v,c(t)v) = f( v ^d t c - [ f C(0), VCG2?([0,T)). 

This type of result is well known; for example, it can be proved along the lines 
of Alt and Luckhaus Q] and Otto [55] (see also [55] and [HJ Theoreme II.5.11]). 

The following lemma is a technical tool adapted to the weak formulation of 
Definition O 



Lemma 2.4. Let 0, be a Lipschitz domain. There exists a family of linear operators 
(K £ )e>o from L 2 (0,T,V) into V(R x R d ) such that 

- for all z G L 2 (0, T, V), 1Z e (z) converges to z in L 2 (0, T, V); 

- for all z G L r (Q) n L 2 (0,T,V), TZ £ {z) converges to z in L r (Q). 

Let us stress that the linearity of TZ e (-) is essential for the application of this 
lemma. It is used to regularise Ui <e , so that one can take 1Z e {ui >e ) as test functions 



in (10 1; for example, a priori estimates for weak solutions and uniform bounds 
on their Galerkin approximations will be obtained in this way. In addition, a 
straightforward application of the lemma is the following uniqueness result: 

Theorem 2.5. Assume (6) and Q. Then there exists a unique weak solution 
(ui,u e ,v) to Problem (JlJ , 1 2 1 , ^ . Moreover, if (ui,u e ,v) is another weak solution 
of Problem ([T]) , ^ , (|3| corresponding to the initial function vq G L 2 (Q), then 

for a.e. te(0,T), IK*) - <K*)l|i»(n) < - «o||i«(n). 

In addition, if Q holds then v depends continuously in L r (Q) on vq in L 2 (f2). 

Continuous dependence of the solution on 7 app , e , ^ e can be shown with the 
same technique, using in addition the Cauchy-Schwarz inequality on (0,T) xTj/ 
and the trace inequalities for H 1 functions. 

Proof. Let C G V([0, T)), C > 0. We take ((t)K s {ui-Ui)(t, x) as test function in the 



first equation of ( 10 ), and ((t)TZ £ (u e — u e )(t, x) in the second equation of ( 10 ). Wc 



subtract the resulting equations and apply the chain rule of Lemma |2.3[ using the 



linearity of TZ £ ( ) and the other properties listed in Lemma 2.4 and subsequently 
sending e — > 0, we finally arrive at 



JJ (Mi(x)(V«i- Vui) ■ (V«i- Vfii) 



+ M e (x)(VM e - Vu e ) • (Vu e - Vw e )J C = 0. 

For a.e. i > 0, we let ( converge to the characteristic function of [0,t]. Thanks to 
the monotonicity assumption Q on h, we deduce 

f(v- v) 2 (t) < f(v- v) 2 (t) + f [ ( h(v) - h(v) )(v-v) 
Jn Jn Jo Jn 

< I K-«o) 2 +2L f f {v~v) 2 . 
Jn Jo Jn 

By the Gronwall inequality, the I? continuous dependence property stated in the 
theorem follows. 

Next, if |8]) holds, from the Holder inequality and the evident estimate 
\v - v\ r < (\v\ + H) r ~ 2 \v - v\ 2 (recall r > 2), 
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we infer that \\v — v\\L r (Q) goes to zero as ||«o — vo||l 2 (0) tends to zero. 

Finally, if vq = vq, not only do we have v = v, but also Ui^ e — Ui tC because of the 
strict positivity of Mj and the boundary/normalisation condition in V . □ 

It remains to prove the regularisation result. 

Proof of Lemma \2.4\ For simplicity we consider separately the two basic cases. 

• Pure Dirichlet BC case. 

Extend z by zero for t ^ (0,T). Take a standard family of mollifiers (p £ ) e> o on 
R d+1 supported in the ball of radius e centred at the origin. Introduce the set 
n e := {x en\ dist (x, dD,) < s}. Take 9 e such that 9 £ € V{tt), e = 1 in Q \ 
< 9 e < 1, and || V9 e \\ L oa {n) < Const/e. Define 

K s (z)(t,x) := (p £ {t,x)) * (9 £ {x)z{t,x)). 

By construction, 1Z E maps L}{Q) to C°°(IR x M d ). From standard properties of 
mollifiers and the absolute continuity of the Lebesgue integral, one easily deduces 
that if z € L r (Q), then z e := lZ e (z) converges to z in L r (Q) as e — > 0. Next, 
consider z € L 2 (0 7 T;V). We have z G L 2 (Q), and thus z e — > z in L 2 (Q) as above. 
In particular, (z £ ) £> o is bounded in L 2 (Q). Similarly, p e * (# e Vz) is bounded in 
L 2 (Q) and converges to Vz = VzinL 2 (Q). Since Vz e = p E *(9 e \7z)+p E *(V9 E z), 
it remains to show that p e * (V6* 6 z) converges to zero in L 2 (Q) as e — > 0. By 
standard properties of mollifiers, it is sufficient to prove that V9 e z — > in L 2 (R d ) 
as e — > 0, which follows from an appropriate version of the Poincare inequality. 

Indeed, in the case d£l is Lipschitz regular, we can fix Sq > and cover O £o by 
a finite number of balls (Oi)i<zi (eventually rotating the coordinate axes in each 
ball) such that for all i £ I, for all e < £o the set fi e fl Oj is contained in the strip 
{^1(^2,^3) < xi < ^i(x2,X3) + Ce} for some Lipschitz continuous function vPj 
on M 2 and some C > 0. Hence by the standard Poincare inequality in domains of 
thickness e, we have \\z(t, OlU^n.) < Ce || Vz(i, Olk^nj- Then 

Jo Jn e e Jo Jn e £ Jo Jo E 

and the right-hand side converges to zero as e — ¥ 0, by the absolute continuity of 
the Lebesgue integral. 

• Pure Neumann BC case. 

We use a linear extension operator £ from V into iJ 1 (IR <i ) such that V H L r (f2) is 
mapped into H 1 (R d ) n i7(M d ). Such an operator is constructed in a standard way, 
using a partition of unity, boundary rectification and reflection (see, e.g., Evans 
P?]). We then define 1Z e by the formula TZ E { Z ) — Pe * (£( z ))- 

• The general case: mixed Dirichlet-Neumann BC. 

It suffices to define il £ := {x e fl | dist (x, Yd) < e}, introduce 8 e as in the Dirichlet 
case, introduce £ as in the Neumann case, and take TZ £ (z) = p e * (9 e £(z)). □ 

Remark 2.6. We have seen that the following space appears naturally: 
E := {(u u u e ) | Ui, e - gi, e e L 2 (0,T,V), v :=Ui-u e £ L r (Q)}. 
Introducing its dual E' and the corresponding duality pairing ((•, •)); we have 

<((x, , (v, i>))) = lim (x, K s y) + (£, n e <p) 

whenever the limit exists. 
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Now, using Remark 2.2 and Lemma |2.4| it is not difficult to recast Problem 
(jTJ) , ([2|) , ([3]) into the following formal framework: 

find (ui,u e ) G E such that (d t v, -d t v) e E' and Q,(|2|,([3]) hold in E' , 

namely, for all (tp,ip) G E, 

' T (((d t v, -d t v), (<p,i>))) 

\Mi(x,VUi) ■ V(p - M e (x,Vu e ) ■ W + h(v)(ip - tp) 

'Sitp- S e l/j) = // I app ((p - ip), 



10 Jr N J 
and for all (tp, ip) € E such that d t tp, d t ip € L°°(Q) and ip(T, ■) = = tp(T, •), 

' T (((d t v, -d t v), ( V ,i>))) = - ff vdt(<p-1>)- I vo(-)&(0,-)-m-))- 



In view of Remarks |2.2| and 2.6 we can apply some of the techniques used by 
Alt and Luckhaus [1] to deduce an existence result from the uniform boundedness 
in E of the Galerkin approximations of our problem (cf. |13|). The uniform bound 
in E is obtained using the chain rule of Lemma |2.3| the Gronwall inequality and 
the assumptions (j6j),([7j on the ionic current. The arguments of the existence proof 
will essentially be reproduced in Section [6j therefore we omit the details here. 

In view of the uniqueness and continuous dependence result of Theorem |2.5| and 
its proof, we can end this section by stating a well-posedness result. 

Theorem 2.7. Assume ^ and |7|. There exists one and only one solution to 
Problem Q, ^ . ^ . If in addition ^ holds, then the solution depends continuously 
in the space E on the initial datum in L 2 (Q). 

3. The framework of DDFV schemes 

We make an idealisation of the heart by assuming that it occupies a polyhedral 
domain of M 3 . We discretise the diffusion terms in ([!]) using the implicit Euler 
scheme in time and the so-called Discrete Duality Finite Volume (DDFV) schemes in 
space. The DDFV schemes were introduced for the discretisation of linear diffusion 
problems on 2D unstructured, non-orthogonal meshes by Hermeline [3H[3S] and by 
Domelevo and Omnes . They turned out to be well suited for approximation of 
anisotropic and heterogeneous linear or non-linear diffusion problems. 

Our application requires a 3D analogue of the 2D DDFV schemes. Three versions 
of such 3D DDFV schemes have already been developed; we shall refer to them as 
(A), (B) and (C). We refer to [1H1[M] for version (A); version (B) that we describe 



in Section 3.2 below was developed in [33 and [3 IH [2]; we refer to [5T] for 
version (C). In this paper, we show the convergence of any of these schemes, using 
only general properties of DDFV approximations. 

3.1. Generalities. In the 3D DDFV approach of [HUGO] (version (A)) and in the 
one of [21 Hill], [3H1I37] (version (£>)), the meshes consist of control volumes of two 
kinds, the primal and the dual ones. Version (C) also includes a third mesh. For 
cases (B) and (C), primal volumes and dual volumes form two partitions of SI, up to 
a set of measure zero. In case (A), the primal volumes form a partition of il, and the 
dual volumes cover 17 twice, up to a set of measure zero. Some of the dual and primal 
volumes are considered as "Dirichlct boundary" volumes, while the others are the 
"interior" volumes (this includes the volumes located near the Neumann part Fjv 
of dfl). With each (primal or dual) interior control volume we associate unknown 
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values for Ui,u e ,v; Dirichlet boundary conditions are imposed on the boundary 
volumes. The Neumann boundary conditions will enter the definition of the discrete 
divergence operator near the boundary; it is convenient to take them into account 
by introducing additional unknowns associated with "degenerated primal volumes" 
that are parts of the Neumann boundary T^. 

We consider the space K x of discrete functions on il; a discrete function u T £ R x 
consists of one real value per interior control volume. On K T an appropriate inner 



product 



is introduced, which is a bilinear positive form. 



Both primal and dual volumes define a partition of into diamonds, used to 
represent discrete gradients and other discrete fields on f2. The space (Hi 11 ) 3 of 
discrete fields on O serves to define the fluxes through the boundaries of control 
volumes. A discrete field Ai T £ (M 3 ) 3 on f2 consists of one M 3 value per "interior" 

diamond. On (R^) 3 an appropriate inner product <^-, ■, J is introduced. 

A discrete duality finite volume scheme is determined by the mesh, the discrete 
divergence operator div s \ : (M 33 ) 3 — > W* obtained by the standard finite volume 
discretisation procedure (with values s x given by the Neumann boundary condition 
on IV)) and by the associated discrete gradient operator. More precisely, the 
discrete gradient operator Vjr : M x — > (M 33 ) 3 is defined on the space of discrete 
functions extended by values g % in volumes adjacent to Tr>; it is defined in such a 
way that the discrete duality property holds: 
(11) 



Vv £ M x , VA4 T £ (R s ) , 



M*, V V r 



..dt 

//r> 
< x = on 
Further, 

denotes an appropriately defined product on the Neumann part Tn of 

the boundary and v a% denotes the boundary values on T N of w T . The precise 
definitions of these objects are given below for version (B). 



Here corresponds to the homogeneous Dirichlet boundary conditional t 
Td, and s x denotes the discrete Neumann boundary datum for M T ■ n. 



In 



differ; 



51] and USUI], [331351, the definitions of dual volumes and 
but both methods can be analysed with the same formalism. The construction in 
PT] only differs by its use of three meshes based on three kinds of control volumes. 



. The main difference between the three 



<> 



This also changes the definition of 

frameworks lies in the interpretation of w x £ R T in terms of functions. In each case 
u % £ M x is thought as a piecewise constant function. The three following lifting 
between W x and L 1 (J7) are considered: 
(12) 



f 



for version (A) described in [49j [21] 

for version (B) described in [3HH], [371EH] 

for version (C) described in pTj . 



with v m ° and i>™* representing the discrete solutions on the primal and the dual 
mesh, respectively, and with v m ° (in the scheme of [H]) representing the solution 



on the third mesh. We have for instance v' m ° (x) 
characteristic function of «r), the definitions ofv m 



E 



ll K (x) (with II* the 



are analogous. 



^our notation follows [6]; a slightly different viewpoint was used in [3] [4] [2], where the ho- 
mogeneous Dirichlet boundary data were included into the definition of the space of discrete 
functions defined also on the control volumes adjacent to To = 9Q. 
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In all the three cases, appropriate definitions of the spaces K T , (IR 11 ) 3 , the scalar 
products •,• >| - , - !r j ((•,■) , and the operators div T , V T , lead to the discrete 



duality property ( 11 ) 



3.2. A description of version (B). In this subsection we describe the objects 
and the associated discrete gradient and divergence operators for version (B) of the 
scheme. More details and generalisations can be found in [2]. 

3.2.1. Construction of "double" meshes. 

• A partition of fHs a finite set of disjoint open polyhedral subsets of SI such that 
fl is contained in their union, up to a set of zero three-dimensional measure. 

A "double" finite volume mesh of SI is a triple T = (971°, 97T , 2?) described in 
what follows. 

• First, let 97f^ be a partition of SI into open polyhedral with triangular or quad- 
rangular faces. We assume them convex. Assume that <9S1 is the disjoint union of 
polygonal parts To (for the sake of being definite, we assume it to be closed) and 
r^r (that we therefore assume to be open). Then we require that each face of the 
polyhedra in 971^ either lies inside S7, or it lies in Tp, or it lies in (up to a set of 
zero two-dimensional measure). Each K G 971^ is called a primal control volume and 
is supplied with an arbitrarily chosen centre x K \ for simplicity, we assume x K G k. 

Further, we call fXftp (respectively, 5977°) the set of all faces of control volumes 
that are included in Tjv (resp., in Tjj). These faces are considered as degenerate 
control volumes; those of 3971° are called boundary primal volumes. For k G 971^ 
or k G <997i°, we choose a centre x K G k . 

Finally, we denote SUt° := 9Ji^ U 971^ ; 971° is the set of interior primal volumes] 
and we denote by 97? the union 97T U dWT = (97^ U 971^) U <997T\ 

• We call neighbours of k, all control volumes l G 971° such that k and l have 
a common face (by convention, a degenerate volume k G 9J$r N or K G 9971° has 
a unique face, which coincides with the degenerate volume itself). The set of all 
neighbours of k is denoted by JV(if). Note that if l G at(k), then k G a/"(l); in this 
case we simply say that k and l are (a couple of) neighbours. If k,l are neighbours, 
we denote by k|l the interface (face) 8k n 8l between k and l. 

• We call vertex (of 97^) any vertex of any control volume k G 9JI^ . A generic 
vertex of 97^ is denoted by x K * ; it will be associated later with a unique dual 
control volume K" G 971*. Each face k\l is supplied with a face centre x KlL which 
should lie in k|l (the more general situation is described in [2]). For two neighbour 
vertices x K * and x L * (i.e., vertices of 97i° joined by an edge of some interface k)z, or 
boundary face), we denote by x K * ]L * the middle-point of the segment [x K *,x L *]. 

• Now if k G 97t° and l G M{k), assume x K -*,x L * are two neighbour vertices 
of the interface We denote by r£» Jf . |i , the tetrahedra formed by the points 
x K , x K * , x KlL ,x K * lL * . A generic tetrahedra t£, k *„, is called an element of the mesh 
and is denoted by t (see Figure [T]); the set of all elements is denoted by T ■ 

• Define the volume k* associated with a vertex x K - of 97i^ as the union of all 
elements t G T having x K * for one of its vertices. The collection 97T of all such k* 
forms another partition of Q. If x K * G SI U Tjv, we say that K" is an (interior) dual 
control volume and write k* G 971*; and if x K * G Tp, we say that k* is a boundary 
dual control volume and write if* G <997T\ Thus 97T = 97T U 997T*. Any vertex of 
any dual control volume k* G 971* is called a dual vertex (of 97f ). Note that by 
construction, the set of vertices coincides with the set of dual centres x K * ; the set 
of dual vertices consists of centres x K , face centres x KlL and edge centres (middle 
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Figure 1. One of the twelve elements in diamond d k{l with 
triangular base ifjL 



points) x K * ]L *. Picturing dual volumes in 3D is a hard task; cf. for version (A) 
and [21] for version (C). 

• We denote by A/"* (if*) the set of (dual) neighbours of a dual control volume if*, 
and by if*|L*, the (dual) interface Ok* fl dh" between dual neighbours if* and l*. 

• Finally, we introduce the partitions of SI into diamonds and subdiamonds. If 
if, l £ are neighbours, let H K be the convex hull of x K and if|L and H L be the 
convex hull of x L and if|L. Then the union H K U H L is called a diamond and is 
denoted by d k{l . 

If if, l £ 9JI° are neighbours, and x K *,x L * are neighbour vertices of the cor- 
responding interface if|L, then the union of the four elements t£* k *, l *, t*, K * lL *, 
t k* k*\l*i anc ^ t l* k*\l* * s called subdiamond and denoted by S^ 1 .^,. In this way, 
each diamond D KlL gives rise to I subdiamonds (where / is the number of vertices 
of ifjL) ; cf . the next item and Fig. [2] Each subdiamond is associated with a unique 
interface ifjL, and thus with a unique diamond d k]l . We will write sCDto signify 
that s is associated with d. 

We denote by D,& the sets of all diamonds and the set of all subdiamonds, 
respectively. Generic elements of £),€5 are denoted by d,s, respectively. Notice 
that 2) is a partition of a subdomain of SI (only a small neighbourhood of Tn in fl 
is not covered by diamonds). 

• (See Figure [2]) The following notations are only needed for an explicit expression 
of the discrete divergence operator (and also for the proof of the discrete duality 
given in [2] ) . It is convenient to orient the of each diamond d. Whenever 
the orientation is of importance, the primal vertices defining the diamond will 
be denoted by x KQ , x K(S in such a way that the vector x KQ x K ^ has the positive 
orientation. The oriented diamond is then denoted by d k&]k ® . We denote by S K&K(S 
the corresponding unit vector, and by (i K0jJCffi , the length of x Kq x k ^ . We denote by 
n Ke \K 9 the unit normal vector to if |if® such that n KQ]K(S ■ e K(SjK@ > 0. 

Fixing the normal ft K@lKfS of if Q |if® induces an orientation of the corresponding 
face if Q |if e , which is a convex polygon with I vertices (we only use / = 3 or 4): we 
denote the vertices of k q \k 9 by x,<*, i £ [M], enumerated in the direct sense. By 
convention, we assign x I ^ i := x^. We denote by e R * )J< * i the unit normal vector 
pointing from 1^ towards %j<*) and by <^k*,k* +1 j the length of i^i^j . 
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Figure 2. Primal volumes, diamond; subdiamond and zoom on it 
(Version (B) of 3D DDFV mesh) 



To simplify the notation, we will drop the if's in the subscripts and denote the 
objects introduced above by d 0ie ,7%) iffl and by x* ,q* i+1 ,d* i+1 whenever 

D K e \K 9 j g g xec [_ We also denote by x* i+1 the middle-point Xk*|k* +1 of the segment 
], and by x* , the centre x Kq]k@ of -R" [fs" e . 

• For a diamond d = d k ® ik ® , we denote by Proj D the orthogonal projection of 
R 3 onto the line spanned by the vector e K@ ,K 9 ! we denote by Proj^, the orthogonal 
projection of R 3 onto the plane containing the interface K^pc®. 

• We denote by Vol(A) the three-dimensional Lebesgue measure of A which can 
stand for a control volume, a dual control volume, or a diamond. In particular, for 
k e 9ttr N , Vo1(k) — 0: these volumes are degenerate. For a subdiamond s = S 1 ^^? , 

y y ' 

we have the formula VbJ(s) = ^ ( x e xl , x*^ x*^ ,x*x* +1 ). Note the mixed product 
is positive, thanks to our conventions on the orientation in d KqIK<s and because we 
have assumed that x* £ -K©^®- 

Remark 3.1. Diamonds permit to define the discrete gradient operator, while 



subdiamonds permit to give formulas for the discrete divergence operator (see ( 13 ) 



( |14| ) and ([17]), ([21]) below, respectively). 

In the context of 2D "double" schemes, introducing diamonds is quite standard 
(see, e.g., [6 ( 26 ). Subdiamonds are "hidden" in the 2D construction : they actually 
coincide with diamonds. 



3.2.2. Discrete functions, fields, and boundary data. 

• A discrete function w' s on fl is a pair (w m ° , w m *) consisting of two sets of real 
values w m ° = {w K ) Keon o and w m * = (w K * ) K *<z 9n * . The set of all such functions is 
denoted by E T . 

• A discrete field M"* on fHs a set (IFn) of vectors of R d . The set of all discrete 

fields is denoted by (R d )" . If M' z is a discrete field on f2, we assign M s — A4 D 
whenever s C d. 

• A discrete Dirichlet datum j x ouTb is a pair (g amt ° , g ™*) consisting of two sets 
of real values g om ° = (g K ) K £ B mt° and .g ™* = (g K . ) K . es£0t » . In practise, g K (resp., 
g^' ) can be obtained by averaging the "continuous" Dirichlet datum g over the 
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boundary volume K C Trj (resp., over the part of Tjj adjacent to the boundary 
dual volume k*); if g is continuous, the mean value can be replaced by the value of 
g at x K (resp., at x K *). We refer to [B] for details. 

• A discrete Neumann datum s x on Tjv is a set of real values (s K ) Kem o . 

In practice, s K can be obtained by averaging the "continuous" Neumann da- 
tum s over the degenerate volume if c Tjv- In the case Tr> = 0, one should be 
careful while using approximate quadratures to produce s x from s. Indeed, some 
compatibility conditions between Neumann data and source terms may arise while 
discretising elliptic equations (this is the case of system ([I]), because the difference 
of the two equations of the system is an elliptic equation, and the compatibility 
condition ^ is needed for the solvability of the system). The compatibility condi- 
tion, expressed in terms of J r s i e , should be preserved at the discrete level. This 
is the case for the above choice of s % : indeed, we have the equality J r s — J r s x . 

3.2.3. The discrete gradient operator. 

• On the set M x of discrete functions iu x on f2, we define the discrete gradient 
operator Vjif-] with Dirichlet data g T on Trj: 

(13) V x * : w* £ R T h+ V g W - (V^^^ G (R d r, 
where the entry VbW 1 of the discrete field V x ui x relative to d 

_ i Proj D (V D w; x =— %,, ffl , 

(14) V D ro T is s.t. < rf ,® 



where 



Proj;(V DW T )= VF(-), 

F(-) is the affine function from K 3 to K that is constant in the direction 
Hq >ffl orthogonal to K|i and that is the ad hoc affine interpolation (namely, 



( 15 ) below) of the values w* at the vertices x*, i = 1, . . . , I, of 
• for the vertices of d lying in fl U Tjv, w Q =w K@ , w 9 = w K@ , w* — Wk*, etc. 
(we use the simplified notation in the diamond d = d k&]K(s , as depicted in 
Figure [2]). For the vertices of d that lie in r^, the values of g" 1 are used: 
e.g., if x Kq £ Tjj, then we set w Q := g KQ in the above formula. 

Clearly, if I = 3 there is a unique consistent interpolation of the values w\, 11)2, W3. 
For I > 4, no consistent interpolation exists, and we choose the linear form in w* 
that leads to the expression 
(15) 

2 1 > 

Proj;(V D w T ) = — — ==f-===f- ^(v&l - W?) [%,e X^]; 

Z^i=l \ %,e i X 0& X ij+l ' X i X i+1 I *=1 

it is shown in [5] that this choice is exact on affine functions, and that it leads 
to the discrete duality formula. For an explicit formula of V D iu x , note that p — 
Proj D ( VuU> x ), p* = Proj^( V D u> x ) are given; then one expresses V D iu x as 
(16) 

w * 1 f VoJ ( s <i<* ) , w , 1/ * * N r — ► ~ — rM 



Remark 3.2. In (14 1, the primal mesh serves to reconstruct one component 



of the gradient, which is the one in the direction . ffi . The dual mesh Wt serves 



to reconstruct, with the help of the formula ( 15 ), the two other components, which 
are those lying in the plane containing if |if ffi . The same happens for version (A) 
of the scheme. On the contrary, version (C) only reconstructs one direction of the 
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discrete gradient on the mesh OJt* , while the third direction is reconstructed on a 
third mesh that we denote by 

Remark 3.3. We stress that our gradient approximation is consistent (see [2] for 
the proof) . Indeed, let w @ , w 9 , (w* i+1 )\ =1 be the values at the points x Q , x 9 , i )'=i , 
respectively, of an affine on d = d Kq]K(b function w. Then V D w x coincides with 
the value of Vw on d. 

3.2.4. The discrete divergence operator. 

• On the set (R^) 33 of discrete fields we define the discrete divergence operator 
divjr[-] with Neumann data s T on IV: 

(17) 

divjt :M T e {R d y h> divJrM* 

where the entries d\v K M. % (for if g SEHq) and div K *A4 T of the discrete function 
div x .A/f x on f2 are given by 

Vif6^, div K M* = — ^ ^ f M D -n Kl 
(lg) VoJ W De3 :^0^ D 

VreaiT, div K ,^ x = — ^- V / M D -n K ,, 

Vol[K*) * ' /«^*nr> 

where (resp., u K *) denotes the exterior unit normal vector to k (resp., to k"). 
Further, for k G %rfr N > we mean that n K points inside Q, and we adapt the following 
formal definition: 

t ig \ V k G , Voi(if )div K 7W x :=M D - n K + s K 

for the diamond d such that 13 fl IV = -K"; 



thus, although Vbi(if ) is zero, in calculations we only use the products Vbj(ii")diy Rr 



r 



which arc well defined thanks to convention ( 19 ). In practise, the discrete equations 



corresponding to volumes of 9J$r N will always read as 
(20) M D ■ n K + s K = 0, 



Notice that the values of the Neumann data s T only appear in the convention ( 19 ) 
for the degenerate primal volumes K C Tjv; at the same time, in the volumes k* 
adjacent to Tjy the data s x are taken into account indirectly. Namely, let K* be a 
dual volume adjacent to IV, and let d be a diamond intersecting K' and adjacent 
to IV; then the value M.d • n K used for the definition of div^* M * is linked to the 



data s x via equations (20) 



The formulas (18) are standard for divergence discretisation in finite volume 
methods; their interpretation is straightforward, using the Green-Gauss theorem. 
The consistency of the discrete divergence operator (in the weak sense) can be 



inferred by duality from the one of the discrete gradient operator (Remark 3.3) and 



from the discrete duality property (11); see Proposition |5.1[ iii) and [2]. 

For the explicit calculation of the right-hand sides in (|18[), one can further split 
diamonds into subdiamonds. In a generic subdiamond, we use the following nota- 



tion. Consider s G &; it is associated with a unique oriented diamond which we 

,K"0|K<B 



denote d k ^ k<s , so that s is of the form s = S^S? . In order to cope with the vector 
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IS 



iven s — 




|0, if if 


= K @ 


11, if k 


= K(B 



we define 



0, if if* — if 

1, if if* = if 



i+1 



For k £ 9Jt°, we denote by v(if) the set of all subdiamonds s € © such that 
k fl s / 0. In the same way, for if* £ 9Ji* we define the set v*(if*) of the 
subdiamonds intersecting if*. Then, using the notation (•,-,•) for the mixed product 



on R , we can express formulae (18) as 



Va'<e9J& div^A^ = 



1 



(21) 



Vif-emr div A ..x T = 



2Vo1(k) ^ ^ S ( "^ s ' ^SffiV' X « ^i+l' 

° W sEv(if) 
1 



2Voj(if) 



^ (-l) e s (TWs,^,^^). 



SGv*(jf*) 



In (21), each subdiamond s in v(if) (or in v*(if*)) has the form s 



with some if , K @ , if* , if * +1 ; the notations , e^'* , x e , % , x* a , , x* , refer to 



,-R<sl«e 



(see Figure 2 ) . Details can be found in [2] . 



Remark 3.4. In practise it is not necessary to calculate the discrete divergence; 
indeed, with the help of the duality property, one can express the discrete system 
of equations in the dual form, where the calculation of the discrete divergence of 
the solution is replaced by the calculation of the discrete gradient of a test function. 
Thus, as for the so-called mimetic finite difference methods, in order to formulate 
a DDFV scheme it is sufficient to calculate discrete gradients. This amounts to the 



"discrete weak formulation" (|27|) we use in the sequel 
3.2.5. The scalar products 



and discrete duality. 
Recall that R' r is the space of all discrete functions on fi. For w' T , u T € R T , set 

= 3 X! Vo K k ) w i< v k + j ^ VoJ(k*) w k ,v k *. 

k£9ji° k* £ an* 

Recall that (R 3 )® is the space of discrete fields on D,. For M % , Q % € (R 3 )®, set 
[M x , G T } = Vo1(d)Mo ■ Go- 



• Recall that in (12), for version (B) of the scheme, given a discrete function w T , 
we set 

Then the function v d ~* £ L°°(r N ) can be defined as the trace of d 1 on Fjv- This 
means, v 9% (x) := -^v K + ^v K * where for H 2 -a.e x £ Tjv, if and if* are uniquely 
defined by the fact that x £ if D if* . 



Finally, for 



we simply use the L? scalar product on T 



N ■ 



Now a straightforward adaptation of the proof of the discrete duality property 
in [5] yields the desired discrete duality property (11). 
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4. The DDFV schemes and convergence results 



The time- implicit DDFV finite volume schemes for Problem ([I]) , ([2| , ^ can be 
formally (up to convention (19)) written under the following general form: 



(22) 



find ((uf' n , u*> n ,v*' n )) C (i? x ) 3 satisfying the equations 

V / n=l,....N 



,x,n+l 



- V % < 



At 



-div|,„ +1 [Mf V£,n+lU 



x,n+li 



LX,n+l _ rx,n+l _ n 



,x,n+l 



- V % ' 



At 



div %, n+1 [M* Vl, n+1 u^ n+1 ] + h 



x,n+l 7-x,n+l r\ 

app u ' 



,x,n+l 



x,n+l 



,x,n+l 



(23) 



)=0, 



For two rigorous interpretations of (22 1, see Definition 4.2 below. 



We normalise u x ' ra+1 by requiring, for all n = 1, . . . , N, 
if r_c = 0, then V" Vol(K)u e>K = 0, 

(24) 



voi(K')u e ^ K , =o, voj(/f°)u eiA .o =o 



(the last condition is only meaningful for the meshing (C)). 

The triple (u*' n+ , uj ,n+1 , i> x,n+1 ) constitutes the unknown discrete functions 
at time level n; Vq and / X pp +1 stand for the projections of the initial datum vq and 

the source term / app on the space of discrete functions. Similarly, 3^' e ™ +1 ,s^'™ +1 
are suitable projections of the Dirichlet and Neumann data gi te ,Si >e , respectively. 
Notice that the boundary data are taken into account in the definition of the discrete 
operators V^x, div s x. The matrices M x e (-) are the projections of Mi, e (-) on the 
diamond mesh. We will mainly work with the mean- value projections; e.g., the 

projection P x on T of Vq would be the discrete function with the entry y ^ K ^ / «o 

J K 

corresponding to a control volume k. For regular functions, the centre-valued 
projection P x can be considered, where the entry v${x K ) corresponds to a volume 



k. We refer to Sections 3.2.2 5.2 for details on the projection operators in use. 

A relation that links h %,n+1 to v T ' n+l closes the scheme; we consider the following 
two choices: the fully implicit scheme, 

(25) /i x <" +1 = P x / l (i/ x ' n+1 (-)), 
and the linearised implicit scheme 

(26) h T ' n+1 = P x ((6(i/ x -"(-)) - L)v T > n+1 {-) -I). 

where P x is the projection operator acting from L 1 ^) into the space of the corre- 
sponding discrete functions; further, u x,n+1 (-) define the piecewise constant func- 
tions reconstructed according to (12 1 from the values v T,n+1 = (v wl °' n+1 , y 3n *> n + 1 \ 
(for versions (A) and (£)) or v T ' n+1 = (v m °' n+1 , v m "' n+1 , v ' m °- n +^) (for version 
(C)). The same convention applies to w x '"(-). We refer to Section 5.7 for a detailed 
description of such discretisation of the ionic current term. 

Remark 4.1. In the discretisation of the ionic current term h(v), the choice 
(25), (26) is made to reconstruct the L 1 function v' x (-) and then to re-project it 
on the mesh X. This is tricky and it may seem unnatural. But we explain in Sec- 



tion 5.7 that this is the way to ensure that the structure of the reaction terms in the 
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discrete equations yields exactly the same a priori estimates as for the continuous 
problem. 

The seemingly simpler choice 



h(v' z,n+1 ) (instead of (25 1) does not have 



good structure properties; we can justify the convergence of the associated scheme 
by adding a penalisation term (cf. [1]) whose role is to make small the differences 



,,71+1 



„n+l 



for K PI K* ^ 0, in the left-hand side of the scheme (22 ) 



Definition 4.2. A discrete solution is a set ( (w^' n+1 , ,n+1 , u T,n+1 



)) (in 

the sequel, we denote it by (uf' &t , mJ' a *, v'* ,&t )) satisfying the initial data (23), the 
normalisation equations (24), and the closure relation (25) or (26); moreover, it 



should solve system (22) in the following sense: 



■ equalities in (|22j) hold component per component for all entries corresponding to 
primal volumes k £ WX^ and those corresponding to the dual volumes k" £ WT ; 



for the entries corresponding to k € 97^ , convention ( 19 ) is used, that is, the 



equations take the form w x >™+ 1 — (uf' n+ — ujf ,,l+1 ) = and 

(M ite ) D V D M T ' n+1 ■ n K + (s^ e )^ +1 =0 for flffD such thatD<T\T N = k; 



Equivalently, (uf' At , ,At , v ' r ' At ) is a discrete solution ifv % ' At = 
and for all (p T € R T , for all n G [0, N] the following identities hold: 



x,At x At 
u A ' — uf' 



(27) 



1 

At 

1 IT 
At 



v x,n+l_ v x,n )(p % 



x,n+l 



v*> n ,<p* 



+ fM?V^ +1 



x,n+l 



- I. 



x,n+l 
app 3 



= o, 



x,n+l 



O 

x,„x 



> v > 



x,n+l 



-7, 



x,n+l 
app ) 



= 0. 



Notice that the equivalence of the two above formulations of (22) is easy to 



establish; namely, the discrete duality property ( 11 ) is used together with the choice 



of discrete test functions i/j x that only contain one non-zero entry. 



The existence of solutions to the discrete equations is obtained in a standard way 
from the Brouwer fixed-point theorem and the coercivity enjoyed by our schemes; 
the uniqueness proof mimics the one of Theorem |2.5| More precisely, we have 



Proposition 4.3. Assume ([6]),([7|. Whenever At < jj-, for all given boundary 
data satisfying ^ (ifTrj = 0, we add (24 1 to the scheme) and for all given initial 
data (23) there exists one and only one discrete solution to the scheme (22), (25 1; 



likewise, there exists one and only one discrete solution to the scheme (22), (26). 
Moreover, for fixed boundary data, the discrete L 2 contraction property holds for the 



v %,At component of the solution of the fully implicit scheme (22),(23),(25): Indeed 
for alln€[0,N], 



(28) 



v x,n+l_Q%,n+l v x,n+l_ -x,n+l 



< L(«+l)At 



w x,0_-x,0 x,0_~x,0 



Remark 4.4. Let us point out that the fully implicit scheme leads, at each time 
level, to a nonlinear system of equations, and to compute the solution given by 
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Proposition 4.3 (or, rather, a reasonable approximation to it) we can use the fol- 
lowing variational formulation of the scheme: 



at the time level n, minimise over 



the functional 



J[u? 



,«f := 



1 

2 At 



1 

A/ 



Mf W»+i«*, V*»,»+itt a 



1,71 + 1 



9c 



,x,n+l 



,9s 



rs,n+l 
app 



where v' 1 := uf — u'* , and H : z i— » /j(s) ds is the primitive of /i 



(in the case = 0, the constraint p4| should be added on the domain of the 
functional J). Similarly to the argument in [B], it is checked from the discrete 
duality formula and from formula (40) in Section 5.7 that the scheme (22 1 is the 



Euler-Lagrange equation for the above problem. From the properties of h(-) it 
follows that for At < Jr, we are facing a minimisation problem for the convex 
coercive functional J. Thus descent iterative methods can be used for solving the 



discrete system ( 22 ) at each time step 



Now we can state the main result of this paper. 

Theorem 4.5. Assume ([6|,([7| hold with some r > 2. Assume that the family 
of meshes satisfies the regularity assumptions (29 1,(30 ), (31 \ (and the analogous 
restrictions on the mesh 9Ji°, for version (C)) stated in Section\5.l\ Then 



(i) the sequence of solution s (u f' *(•), ujf' At (-), v' T ' At (-)) to the fully implicit 
scheme (22), (23 1, (25), (12) converges, as the approximation parameters 
Ax, At tend to zero, to the unique solution (itj, u e , v) of Problem (jTJ) . ([2]) , (J3j) ," 
the convergence is strong in L 2 (Q) x L 2 (Q) x L r {Q). Moreover, the discrete 

gradients converge to (Vui, Vit e , Vw) strongly in (L 2 (Q)) 3 ; 

(ii) For any r < 16/3, the statement analogous to (i) holds for the discrete 
solutions of the linearised implicit scheme (22), (23), (26), (12). 



If Yd = 0, the constraint ( 24 ) should be added to the equations of the scheme. 



In the same vein, the standard 2D DDFV construction can be applied to problem 
(jTJ) , ([2]) , ([3]) on 2D polygonal domains. The convergence result of Theorem 4.5 i) 
remains true, and the one of Theorem |4.5[ ii) extends to all r < 6. We stress that 
the realistic case r = 4 is covered by our convergence results. 



5. Discrete functional analysis tools for DDFV schemes 
For a given mesh X of fl as described in Section |3j the size of X is defined as 

size(X) := max < max L diam(i<") , max diam(i<'*) , maxdiam(_o) > . 

If the assumption x K G k is relaxed, diam(if) must be replaced with diam(Kri{x Ji: }) 
in the above expression. 

In what follows, we will always think of a family of meshes such that size(T) 
goes to zero. 
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5.1. Regularity assumptions on the meshes. 

In different finite volume methods, one always needs some qualitative restrictions 
on the mesh X (such as, e.g., x K £ k, or the convexity of volumes and/or diamonds, 
or the mesh orthogonality, or the Delaunay condition on a simplicial mesh). For 
the convergence analysis with respect to families of such meshes, it is convenient 
(though not always necessary) to impose shape regularity assumptions. These 
assumptions are quantitative: this means that the "distortion" of certain objects 
in a mesh is measured with the help of a regularity constant reg(T), which is finite 
for each individual mesh but may get unbounded if an infinite family of meshes is 
considered. For the 3D DDFV meshes presented in this paper, there are two main 
mesh regularity assumptions. First, we require several lower bounds on d KL ,d K * L *: 

V neighbours k, l, diam (k) + diam (l) < ieg(T)d KL ; 

V dual neighbours k*,l* 7 diam (.ft") + diam(z,*) < reg{%)d K * L * ; 

V diamonds d with vertices x K , x L and with 

neighbour dual vertices x K * ,x L *, diam (z>) < reg(T) min{d KL , d K * L * }. 

Further, we need a bound on the inclination of the (primal and dual) interfaces 
with respect to the (dual or primal) edges: 



(29) 



(30) 



V primal neighbour volumes if, l, the angle a KtL between x K xl and the 
plane k\l is separated from and 7r, meaning that reg(T) cosa KiL > 1; 

V neighbour vertices x K * 7 x L » of if|L, the angle a*,*^, between x K *x L * 
and x K * lL *x Kl l is separated from and ir, i.e., reg(X) cos a*, L * > 1. 



Also a uniform bound on the number of neighbours of volumes / diamonds is 
useful: 



(31) 



Each primal volume k has at most reg(X) neighbour primal volumes; 
each dual volume k* has at most reg(T) neighbour dual volumes. 

For version (C) of the scheme, we impose in addition conditions on the third mesh 
9Ji°; moreover, the number of vertices of a diamond is restricted by reg(X). Recall 
that for versions (A) and (B) we assume that all diamond has five (— 2 + 3) or 
six (=2 + 4) vertices, because the faces of the primal volumes are assumed to be 
triangles or quadrilaterals; and the convergence results are shown for the case of 
triangular primal faces. 

For versions (A) and (B), when the number I of vertices of a face if|z, exceeds 
three, the kernel of the linear form used to reconstruct the discrete gradient in d k]l 
is not reduced to a constant at the vertices of d 1 ^ . This is a problem, e.g., for the 
discrete Poincare inequality and for the proof of discrete compactness. In general, 
the situation with I < 4 vertices is not clear; for example, the discrete Poincare 
inequality holds on every individual mesh, but it is not an easy task to prove that 
the embedding constant is uniform, even under rigid proportionality assumptions 
on the meshes. The uniform Cartesian meshes is one case with I = 4 that can be 
treated (see [2]), but they are not suitable for the application we have in mind. 

In this paper, for a certain range of values of the power r in (|6|,([7]), we use 
Sobolev embedding inequalities of the discrete H 1 spaces into L q , q > 1; for these 
results to hold, we may also require 

V primal volumes k and interfaces ^,m mL d KL < reg(T) VbJ(_fs"); 

V dual volumes k and interface _kt*|l*, m K * lL *d K + L * < reg(T) Vo1(k"). 



(32) 



5.2. Consistency of projections and discrete gradients. Here we gather basic 
consistency results for the DDFV discretisations. Heuristically, for a given function 
(p on J7, the projection of <p on a mesh X and subsequent application of the discrete 
gradient V -1 should produce a discrete field sufficiently close (for size(T) small) to 
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Vip. Similarly, for a given field M, the adequate projection on the mesh and the 
application of div x to this projection should yield a discrete function close to div M. 
In this paper, we mainly use the mean-value projections. For scalar functions on 
ft, two projections on K x (which has two components, namely the projections on 
and on 9Jf ) are used: 

pI : ^ ( ( i • ( ^) I ) = : ( p ~* ■ p ^ ) < 

v * ((vWU- (^)) if . eSDr )=:(PfV,PrV); 

in case K is a degenerate volume in 0Jt^ N , VoI(k) is zero and we replace the cor- 
responding entry of P x ip by the mean value -f K <p of ip over the face k C T^. 
Similarly, the Neumann data Si yC will be taken into account through the values 
£s i)e for k€ 

Further, if we are interested in the values of ip on the Dirichlet part of the 
boundary, then we use the projection 

v) Kem ? , (f_ ^ K , eaajP ) =: (P*"V , P 9OT V) • 
\ Jk J K*r\T D / 

In particular, the Dirichlet data gi e will be taken into account in this way. For 

M 3 -valued fields on fl, we use the projection on (M 3 ) 1 * defined by 

P T : M h> ( — ^ [ M 

With each of these discrete functions, we associate piecewise constant functions 
of x on f2, on T/j or on T^, according to the sense of the projection; then we can 
study convergence, e.g., of P X A^ to M. in Lebesgue spaces, as size(T) — > 0. For the 
data vo,/app) M^e, <?j,e> Si,e, we need the consistency of the associated projection 
operators (recall that i>o,/app are projected on the meshes 97^ and SUt*, Mj. e are 
projected on the diamonds, gi^ e are projected on the boundary volumes, and Sj e are 
projected on the degenerate interior primal volumes k G 97l£ ). These consistency 
results can be shown in a straightforward way (see, e.g., [5]); for example, we have 
P x /app — > 4 PP in L 2 (ft), and P 9 ' 1 ft , e — > 5l , e in L 2 {Y D ). 

Note that for the study of weak compactness in Sobolev spaces and convergence 
of discrete solutions, the consistency results can be formulated for test functions 
only (and the consistency for div x oP x is formulated in a weak form, except on 
very symmetric meshes). These results are shown under the regularity restrictions 
( 29 ),( 30 ) , (pl]> on the mesh; let us give the precise statements. 

Proposition 5.1. Let % be a 3D DDFV mesh of VL as described in Section^ 
Let reg(T) measure the mesh regularity in the sense of (29 1,(30 ),( 31 1. Then the 
following results hold: 

(i) For all tp G V(U), 

\W- pOT V|| i0 o (n) < C(<p) sizeCX), \\<p- P^IL-pi) < C(ip) sizeCX); 
and for version (C), the analogous estimates hold for ||</> — P'^VH^cwqv 
Analogous estimates hold for the projections ¥™° . 
Similarly, for all M G (f(O)) , 



\M - W % M\\ LOO{ < C(M) size{%). 



(ii) For all ip G T>(Q U T N ), 



V<p- Vj(Pfr)|| ioo(n) < C(^,reg(T))^e(T). 
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(Hi) For versions (A) and (B), assume that each primal interface k\l is a triangle. 
For each M G {V{VL)f and for all w T £ R% , 

¥ % (divM) - div T (P T 7W) , w T < C(7W,reg(X)) size{%) || V T u; x || L i ( n). 

We refer to [5] for a proof of this result. 

5.3. Discrete Poincare, Sobolev inequalities and strong compactness. 

The key fact here is the following remark: 

Assuming (for versions (A), (B)) that each face k|l of SP1° is a triangle, 
one gets the same embedding results on the 3D DDFV meshes (A), (B), (C) 
as the results known for the two-point discrete gradients on 2JT and on WT . 

Indeed, for variants (A), (B) it has been already observed in the proof of Proposi- 
tion 5.1 hi) that the restriction i = 3 on the number I of dual vertices of a diamond 



d k & \k 9 a vj ows f or a control by | V D w x | of the divided differences: 

\w m — wj\ i„ t i \w* u —w*\ , 
(33) ^ < | V D w T |, i-^ ^ < | V D w*\ 

(here i = 1,2, 3 and by our convention, w\ := w*, d^^ := di^; cf. Figure [2|. For 
version (C), this kind of control is always true for the divided differences along 
the edges of any of the three meshes. Consequently, for a proof of the different 
embeddings, we can treat the primal and the dual meshes in X separately, as if our 
scheme was one with the two-point gradient reconstruction. 

First we give discrete DDFV versions of the embeddings of the discrete Wq ,p {£1) 
spaces, where we refer to the embedding into L p (fl) (the Poincare inequality), into 
LP* (O) with p* := p < 3 (the critical Sobolev embedding), as well as the 

compact embeddings into L q (Q) for all q < +oo or q < p* . 

Proposition 5.2. Let X be a 3D DDFV mesh of fi as described in Section [5| 



Let reg(X) measure the mesh regularity in the sense (30) and (32). Assume (for 
versions (A) and (B) ) that each primal interface k\l is a triangle. 
Let w T £ Rjf . Then 

\\w m "\\mn),\\w m *\\ LH n) < C(n,reg(SO)||V*w*|U» ( n}. 

Moreover, 

K m lL« ( n), \\iv™*\\Le(n) < C(fi,reg(X)) || 



Notice that for the Poincare inequality (the first statement), assumption (32) is 



not needed, cf. [7] for a proof. Actually, with the hint of Lemma 2.6] the Sobolev 



embeddings for g < 2 x 1* = 3 can be obtained without using (32) 



The statement follows in a very direct way from the proofs given in [551 E01 US] • 



Because of (33), the assumption that the primal mesh faces are triangles (i.e., 
I = 3) is a key assumption for the proof. In some of the proofs in these papers one 
refers to admissibility assumptions on the mesh (such as the mesh orthogonality 
and assumptions of the kind tl \x K — x L \ < reg(X)|:n K — x K]L \" , see (28j[29]), yet, as 
in [B] (where the proof of the Poincare inequality is given for the 2D case), these 
assumptions are easily replaced by the bounds 

m KlL d KL < C*(reg(X)) min{ Vo1(d^ l ), Vo1(k), Vo1(l)}, 

m K * lL *d K * L , < C*(reg(X)) iam{Vol( D ^,), VoI(k'), Vo1(l*)} 



that stem from the mesh regularity assumptions ( |32| ) and ( 30 1 . 

The embeddings of the discrete W 1,p (ft) space contain an additional term in the 
right-hand side, which is usually taken to be cither the mean value of w x on some 
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fixed part T of the boundary <9f2 (used when a non-homogeneous Dirichlct boundary 
condition on T is imposed), or the mean value of iu x on some subdomain w of Q (the 
simplest choice is ui = CI, used for the pure Neumann boundary conditions). Let us 
point out that the strategy of Eymard, Gallouet and Herbin in [25] actually allows 
to obtain Sobolev embeddings for the "Neumann case" as soon as the Poincare 
inequality is obtained. For the proof, one bootstraps the estimate of J n \w' z \ a . 
First obtained from the Poincare inequality with a = 2, it is extended to a = 
2 • 1* = 2 1 = 3 with the discrete variant [251 Lemma 5.2] (where one can exploit 



(34)) of the Nirenberg technique. In the same way, the bound of J Q ju^l" is further 

extended to a = 2(1*) 2 = 2(|) 2 and so on, until one reaches the critical exponent 
2* = 6. The details are given in [5]. Moreover, the Poincare inequality for the 
"Neumann case" (i.e., the embedding into L 2 (fl) of the discrete analogue of the 
space {u £ H^fl) | J u u = 0}) and for the case with control by the mean value on a 
part of the bou ndar y, was shown in |29j , |30j . Thus we can assume that the analogue 



of Proposition 



5.2 



with the additional terms | J n w m ° \, \ J n 



w ' 



| jp^y J Td w m ° | , | -pp^r Jp o w' m * | in the right-hand side of the estimates is justified. 

Notice that the compactness of the sub-critical embeddings is easy to obtain by 
interpolation of the L 6 embedding with the compact L 1 embedding derived from the 
Helly theorem (indeed, the L 1 estimate of V 1 !!) 1 can be seen as the BV estimate 
of the piecewise constant functions w™° and w m *). 

Finally, notice that the same arguments that yield the Poincare inequality with 
a homogeneous boundary condition also yield the trace inequality 

( 35 ) IK OTO L(r N ) ^ C ( r N,n,reg(T)){\\w^\\ L2m + || W|| L2(n)> 

(the inequalities on the mesh Wt and, for the case (C), on the mesh *Xft are 
completely analogous). These inequalities are useful for treating non- homogeneous 
Neumann boundary conditions on a part r^v of dQ. 



5.4. Discrete W 1:P (ft) weak compactness. In relation with Proposition 5.2 ji), 
let us stress that there is no reason that the components w™° h , w Dyt * h of a sequence 
(w Th )h of discrete functions with LP bounded discrete gradients converge to the 
same limit. Counterexamples are constructed starting from two distinct smooth 
functions discretised, one on the primal mesh SW, the other on the dual mesh 9Jt*. 
However, the proposition below shows that in our 3D DDFV framework, we 



w %h — (u!™° h , w' m " h , w m<>h ) coincides with the limit of (12) 



Proposition 5.3. 

(i) Let w Th € R' T '' be discrete functions on a family {%h)h of 3D DDFV meshes of 
fl as described in Section^ parametrised by h > size(T^). Assume ^ and 
let g T = P T <?, for some fixed boundary datum g € _ff 1 (il) . Assume that the family 
(V^W*") h ^ Qfhmam] is bounded in L 2 {fl). 

Assume (for versions (A) and (B) ) that each primal interface k\l is a triangle. 
Assume that sup^Q ^^j reg(T/,.) < +oo, where reg(T/ t ) measures the regularity 



of% h in the sense (29), (301,(31) and (32) 



According to the type of 3D DDFV meshing considered, let us assimilate w %h 



into the piecewise constant functions w Th (-) defined by (12); furthermore, let us as- 
similate the discrete gradient Vj^ h w Xh to the function ( Vj£ h w Th ) (•) on Q. Then 
for any sequence (hi)i converging to zero there exists w € g + V such that, along a 
sub-sequence, 

w Th i (•) converges to w strongly in L 2 (il) (in fact, in L q (fl), q < 6) 
and ( V Thi u' ! ' h i)(-) converges to Vw weakly in L 2 (£l). 



(36) 
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(ii) If To — 0, and if the additional assumption of uniform boundedness of 



1 



m„„ajt° 



Vbj(O) 



1 



171 gji* 



Voi(fi) 



is imposed (with the analogous bound on the mesh 971° for version (C)), then (36) 
holds with w G (f2). 

Let us illustrate the DDFV techniques by giving the ideas of the proof. We 
justify (i) for the case of the meshing described in Section [3] and the homogeneous 
Dirichlet boundary condition 011^ := dfl. The case of a non- homogeneous Dirich- 
let condition is thoroughly treated in 0, for the 2D DDFV schemes. The case of 
Neumann boundary conditions is the simplest one. 



Proof of Proposition \5.3[ The strong compactness claim follows by the compactness 
of the subcritical Sobolev embeddings of Proposition 5.2 The weak L? compactness 
of the family ( W Th w Th ) h is immediate from its L 2 (f2) boundedness. Thus if w is 



the strong L limit of a sequence w %h 



and x is the weak 



L 2 limit of the associated sequence of discrete gradients W Th w' Zh , it only remains 
to show that x = ^ w m the sense of distributions and that w has zero trace on 
dfl. These two statements follow from the identity 



(37) 



w div M = 0, 



X-M 

that we now prove. We exploit the discrete duality and the consistency property of 

and write the discrete duality 



Proposition (5.1 1 (i), (hi). 



Take the projection ¥* h M G (K 3 ) 1 " 1 , 
formula 



(38) 



According to the definition of 



h M 



w ' 



div Xh P' Th .M 



0. 



, the first term in ( 38 ) is precisely the integral 



over f2 of the scalar product of the constant per diamond fields V T ' l w ! ' h and P x,l .M. 
By Proposition (5.1| (i) and the definition of x, this term converges to the first term 
in (37) as h —> 0. Similarly, introducing the projection P Xh (div.M) of d\v A4 on 

R Th , from the definition of • , • , Proposition (5.1 )(i) and the definition of w %h 



in (|12|), we see that, as h - 
P T (divA^) 



w ' 











.1 n 





(lim w m ° h 



(lim w 



div.M 



w div M. . 

It remains to invoke Proposition (5.1)(iii) and the L 1 (i7) bound on V rh u; x ' 1 to 
justify the fact that 



lim 

/i-s-0 



,i h div Tft P' 3rh A^ 



= lim 

SI h^Q 



For a proof of (ii) use the versions of the compact Sobolev embeddings with control 
by the mean value in J7, and use test functions Ai compactly supported in 51. □ 

5.5. Discrete operators, functions and fields on (0, T) x f2. We discretise our 
evolution equations in space using the DDFV operators as described above. In this 
time-dependent framework, analogous consistency properties, Poincare inequality 
and discrete L p (0, T; W 1,p (£l)) compactness properties hold. 
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To be specific, given a DDFV mesh X of ft and a time step At, one considers the 
additional projection operator 

S At : / h+ (/ n ) n6[1 ^ At] C £ X (fi), /"(or) := - / /(t,x) A. 

A,/ ./ (n— l)Ai 

Here / can mean a function in L 1 ((0,T) x ft) or a field in (L 1 ((0,T) x ft)) 3 . The 
greatest integer smaller than or equal to T/eA is denoted by A^ t . 

We define discrete functions w x,At £ (R*)^' on (0,T) x ft as collections of 
discrete functions w' T '" +1 on ft parametrised by n £ [0, A At ] PlN. Discrete functions 
w s,At g ( K Tjiv at Qn ^ 0)T ^ x q and discrete fields A4 x ' At G (r)^' are defined 
similarly. The associated norms are defined in a natural way; e.g., the discrete 
L 2 (0,T;#g(ft)) norm of a discrete function w' T ' n £ (R%) NAt is computed as 

Er>ii vvn+i i^cn) 

To treat space-time dependent test functions and fields as in Proposition |5.1[ one 
replaces the projection operators P x (and its components f m ° ,F m "), ¥ dT and px 



by their compositions with S At . Then the statements and proof of Proposition [5T 
can be extended in a straightforward way. 

Also the statements of Proposition |5.3| extend naturally to the time-dependent 
context; one only has to replace the statement (36) with the weak L 2 (0, T; if 1 (ft)) 
convergence statement: 



IV 



?h,Ath 



converges to w weakly in L 2 ((0,T) x ft) and in L 2 (0, T; L 6 (ft)); 



^ % h w th,Ath converges to Vic weakly in L 2 (ft), 

as size(Xfr,) + At/, — > 0. It is natural that strong compactness on the space-time 
cylinder (0, T) x ft does not follow from a discrete spatial gradient bound alone; one 
also needs some control of time oscillations. It is also well known that this control 
can be a very weak one (cf., e.g., the well-known Aubin-Lions and Simon lemmas). 
In the next section, we give the discrete version of one of these results. 

5.6. Strong compactness in i 1 ((0, T) x ft). Below we state a result that fuses a 
basic space translates estimate (for the "compactness in space" ) with the Kruzhkov 
L 1 time compactness lemma (see |41j). Actually, the Kruzhkov lemma is, by 
essence, a local compactness result. For the sake of simplicity, we state the version 
suitable for discrete functions that are zero on the boundary; the corresponding 
Lj oc ([0, T] x ft) version can be shown with the same arguments (cf. [5]), and this 
local version can be used for all boundary conditions. 

Proposition 5.4. Let (u Th ' Ath ) h £ (RQ h ) NAth be a family of discrete functions 
on the cylinder (0, T) x ft corresponding to a family (At^jh of time steps and to a 
family (%h)h of 3D DDFV meshes of ft as described in Section^ we understand 
that h > size{%h) + Ath- Assume that sup feg ( : ^ max | reg(Xfo) < +oo, where reg(Xh) 



measures the regularity of %h in the sense (29) and (30). 

Assume (for versions (A) and (B) ) that all primal interfaces k\l for all meshes 
are triangles. For each h > 0, assume that the discrete functions v Th ' Ath satisfy 
the discrete evolution equations 

for n £ \0,N h ], = div Xh M %h ' n+1 + 

At 

with some initial data v Th '° £ R Th , source terms f^^Ath g (^^h^N &th an( ^ di scre t e 
fields M %h ' Ath £ ((R 3 )* h ) NAt K 
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Assume that there is a constant M such that the following uniform L 1 ((0, T) x f2) 
estimates hold: 



^(lh'~ +1 IL 1( o)+ll- 



, II fm° h ,n+l || , || fin* h ,n+l II 

II 3 llL 1 (f2)" r II Nil (£2) 



lL!(f2) 



< M, 



and 



At || V Th u Th '" +1 . < M. 



Assume that the families (b(u m ° h '°)^ h , ( b(u™* h '°) j are bounded in L 1 (f2). 

Then for any sequence (hi)i converging to zero there exist j3°, (3* £ L 1 ((0, T) x f2) 
sitc/i i/iai, extracting if necessary a sub-sequence, 

b(u m " h i' Ath i) — > 13°, b(u m * h *> Ath *) — > (3*, in ^((O,! 1 ) x fi) as i -> oo. 



Notice that we only use the full strength of Proposition |5.4| to treat the linearised 
implicit scheme. For the fully implicit scheme, more traditional (although not much 
simpler) I? versions of time translation estimates, inspired by the technique of 1 , 
can be used (see [25]). 

5.7. Discretisation of the ionic current term. Consider the general situation 
where w is discretised on a DDFV mesh. Moreover, assume that we also need to 
discretise some scalar function ip{w) (in our context, this is the ionic current term 
h; in general reaction-diffusion systems, ip may represent reaction terms). 

Then we discretise such reaction term on a 3D DDFV mesh of the kind (B) by 
taking, for W — ip(w), 



: = 



(39) 



In other words, 



voi( K n *:*) 



1 2 
3^+3 ^ 

1 



Vo1(k) 



E 



voi( K n «r») 

Voi(ir*) 



u> K and are the mean values 
of the function w x (-) 

With this choice, we have for all w % £ KJ and for all (p^ £ 



^w^ 1 " (■) + |w OT *(-) on k and on if*, respectively. 



For schemes (A) and (C), we use similar projection formulas with the expression 
of u> x (-) given by (12 1; this always leads to the formula 



(40) 



(vnr,^ = / vK(-))^ x (-) 



By the definition (12 1 of w' T (-), the definition of 
we get 



and Jensen's inequality, 



(41) 



Vw x £ 



k*(-)l< 



Finally, notice that such choice of discretisation of the ionic current term does not 
enlarge the stencil of the DDFV scheme used for the discretisation of the diffusion. 
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6. The convergence proofs 



6.1. Convergence of the fully implicit scheme. The proof follows closely the 
existence proof for Problem ([!]) ,([2]), ^ mentioned in Section [5] 

Step 1 (proof of Proposition 14.31 - uniqueness of a discrete solution) . Although 
Remark |4.4| can be used to infer the existence and uniqueness of a discrete solu- 
tion, let us give a proof that contains the essential calculations also utilised in the 
subsequent steps. For the uniqueness and the continuous dependence claim (|28|) we 



reason as in Theorem 2.5 omitting the regularisation step. Namely, using (22 1 for 
two 
get 

(42, e ) 



_ o — o sr j i o ' ■ 

two solutions (^(u^' At , uJ' At , v T ' At )^ and ^(u l x,At , -uJ' At , -0 X ' A *)^ , by subtraction we 



-(-l)^divl,» +1 [M^ e (Vl„ +1 ' 



£,7l+l 



with (-!)*:=!, (-l) e :=-l and/i T >™ +1 = F T h(v T ' n+1 {-)), h T - n+1 = F x h(v 



1,71+1 



(•))• 



For all n, we take the scalar product 



of equations (42 4 e ) with the discrete 



functions ip 



X ._ ( n ,-z,n+l 



— u x,n+1 ), respectively (more precisely, we use the dis- 



crete weak formulations (27 1 with test function <p T ). Then we subtract the relation 
obtained for e from the relation obtained for i. Finally, we use the discrete dual- 
ity property (11) on the divergence terms. Notice that the boundary terms van- 
ish, because both solutions correspond to the same Dirichlet and Neumann data 
" ~ ' A *,sfg At ; in particular, we can use ( [TT| ) because 



9i,\ 



= V,f (Ut n+1 -Ule n+1 ) 



r ,v ,i , • r -j— \ // T,n+1 x,n+l x,n+l ~x,n+l 

further, the terms coming from 1 n are <ls i e — s i e ,u i e —u i e 
The outcome of the calculation is the following equality: 

(43) 



= 0. 



1 

a7 



w i,n+l_^*,n+lJ _ ( v *,n_ fi *,n) j ^x.n+l _ -x,«+l) 



f(Vl,„+lV -Vl,„+l«i' 



, ( V^r.n+i^ - V g x,„ + iW i J 



MJ(V|,„ + 1 U 



X,Tl+l 



x,n+l"\ 



1,71+1 



V^.„+iu, 



x,n+l\ 



,1,71+1 ux,n+l 



).(« 



x,n+l_^x,n+l , i 



0. 



<> 



Then we sum over n £ [0,k], k < N. Using the convexity inequality a(a — b) > 
\{a 2 — b 2 ), the positivity of M x e and the definition of b, using (40) we get 



(44) 



^X,fc+l_£X,fc+^ 



n=0 



< 



„x,0 a*,0\ 



( r " r ~ "' ) , (v % '" — v" 



At 



,1,71+1 



(•)-^' n+1 0l 



Then the second term is non negative, and we get (28) from (41) and the discrete 
Gronwall inequality. 

In particular, it follows that for fixed initial and boundary data there is unique- 
ness of 7j T < fe+1 for all k. Then, returning to (43 ), we find out that there is uniqueness 
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for Vjr,„+iW^g l+1 for all n. We conclude the uniqueness of uj^ +1 using the discrete 

Poincare inequality (and, in the case To = 0, using condition (p4|). 
Step 2 (proof of Proposition 



4.3 



existence of a discrete solution). Regarding 
the question of existence, we reason by induction in n. Using the discrete weak for- 
mulation (|27[) with the test function <p x = uf e , subtracting the equations obtained 



n+1 , V|,„ +1 < 



for subscripts i and e, we find 

l it 



At 



<> 



Mf V|,„ +1 U, 



n 

x,n+l 



Q 
,n+l 



l it 

A* 



app 







■ n + j 




x,n+l 









x.n+1 T,n+1 



Then using the condition > L, property (41 1, the Cauchy-Schwarz inequality 



and the equivalence of all norms on K x , we deduce the a priori estimate 



7 f J t-tx x,n+l t-7 x x,rs+l 



,n+l 



x,n+l 



1 

2~Ai 



- L 



v t,n+l v %,n+l 



< C*(v T < 



x,n+l *,n+l x,n+l 



' app 5 ° 



z,e 



,/,7,f7,r/j). 



The left-hand side allows us to bound the discrete solutions a priori, if Trj ^ 0. 
The case of pure Neumann BC is slightly more delicate. 

We take advantage of the above estimate to apply the Leray-Schauder topological 
degree theorem. Let us look at the most delicate case Td = 0. 

For 9 £ [0,1], we consider the initial data Ovq, the Dirichlet and Neumann 
boundary data 9gi, e and 9si. e , and the source term 9I app . We consider a family T e 
of maps on the space 



Sp:={((u?' n ,u*'», «*'")) C(R*) 3 

L \ / n— 1.....JV 



(24| holds for all n 



x,n+l 



introduce the 



defined as follows. First, given an element in Sp denoted U 
following notation for the expressions in the left-hand side of equations ( 22 ) with 
data scaled by 9: 



al(U 
a e e (U^ At ) : 



X,At\ ~ ^t|•'-'" + l -^> T ■^ 



A/ 



div^x,„ +1 [Mf V-x,„ +1 ^'" +1 ] + h^+i - 9I*£+\ 



At 



6si 



'94 

- div^ ? ,„ + i [Mf V e x g x,„+i<' n+1 ] + /i x ' n+1 - 



app 

rx,n- 
app 



recall that for the volumes k £ ^fXtp N , the convention ( 19 1 applies, so that the 
entry of a® e (U % ' n+1 ) corresponding to the volumes k £ £01^ should be taken 
equal to (M ie )™ +1 • v K + (s i e )™ +1 , where the diamond d is the one with k C 3d. 
Then we define F 6 as the element of K' T xl s x K' T given by (of, of - a*, 7 e ) . This 
definition implies that JF e maps into itself, thanks to the definition of the discrete 
divergence (which ensures the consistency of the fluxes) and to the constraint ([5| 
that is preserved at the discrete level. 

With this definition, it is evident that the zeros of T are solutions of the scheme 



(22 1 with data scaled by 9. The estimate that we have just deduced is uniform in 



9; it provides a uniform in 9 bound on some norm of possible zeros of T in Sp. 
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Therefore from the Leray-Schauder theorem and the existence of a trivial zero of 
J 70 we infer existence of a zero for T , in particular for 8=1. 

Step 3 (Estimates of the discrete solution). We make the same calculation as 



in Step 1, but with uf' e t set to zero; and we sum over n = 1, . . . , k. Using the 
convexity inequality a(a—b) > \{a 2 — b 2 ), the positivity of M.f e and the definition 
of b, using (40 1 we get the identity 

+ 



(fc+l)At 



'7 



(fc+l)At 



(|(V^ At ^ A4 )(0| 2 + |(V|,A t ^' At )(-)| 2 ) 



< 



(fe+l)At 



(fe+l)At 



+ / / (L|^' At (-)| 2 + ^ T ' At (-)) 



4 T P At (-)^' At (-) 



9«.A* ( . ) + a fl«,A* ( . )|tt ««,A* (0 



Using the Cauchy-Schwarz inequality, property (41 1, the trace inequality (351 for 
each components of the solution, and the discrete Gronwall inequality, we deduce 
the following uniform bounds: 

(45) ||« t ' A *(-)IIl~(o,t ; l^)) < C; 

(46) \\v %At {-)\\ L r {Q) <C; 

(47) ll<; At (-)llL 2( Q) + ||(v 9 \ At <; At )(-)|| i2(Q) < C, 

where C depends on reg(T), 7, a, L, I, ||^olU 2 (n), ||4p P ||l2(q), \\9i, e \\L^(o,T;H^a)), and 
ll s i,e||L 2 ((o.T)xr JV )- As usual, in order to obtain (47), the case To — is treated 
separately, using the normalisation property (24) for it e , the L 2 (Q) bound on v and 
the fact that u e = itj — v. 

Step 4 (Continuous weak formulation for the discrete solutions). We take a test 
function 99 6 T>((0,T] x (f2urjv)) and discretise it as follows: 

^i,Ai ._ p* §At^. Qn tlle Di r i c hlet boundary, we take := 0. 



Then we use the discrete weak formulation (27) with test function Atip' T ' n+1 at time 
level n, and sum over n. What we get is 



N 



N 



(48) J2 \\v T ' n+1 -v T ' n , <p % < n+1 



+ Y l At\\h*' n+1 ,<p 



x,n+i 



n=0 



N 



n=0 
N 



T,n+1 



T,n+1 



n=0 n=0 
,x,A* 



The equation for the components it J' is analogous. 



We use summation by parts on the first term in ( 48 1 , the Lipschitz continuity of 
dt<p, the definition of i> x, ° and Proposition 5.1 'i), to see that this term equals 

- [[ v*> At (-)d t <p- I n ip(0,-)+rl( S \zeCZ),At)(l + \\v \\v m ), 
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where r v denotes a generic remainder term such that r v (size(X), At) — > as 
size(X), At — > 0. Thanks to (40), the second term in (48) is merely 



N 

J2*t I h(v*' n+1 (.))<p*> n+1 (.)= II h(v*> At (.))<p*> At (.) 

%,Atf \\ ,„ , „ („\„ n r<T\ A .*\ r<<^ n -iT\\\\'ui„,^^\ 



^ h(v*>**(-)) ^ + 7>(size(X), At) C(reg(<Z))\\h(v*' At )\\ L1{Qy 
Because the discrete gradients are constant per diamond, thanks to th e de finition 



of M^™ -1 " 1 and to the discrete gradient consistency result of Proposition 5.1 4i), the 
third term in ( 48 ) is equal to 



N 

E 

n=0 



At ( M^V|,„ +I ^' n+1 (-)- V ^'" +1 (-)= // AtU f' At (-)- 



, (size(X) , At) C7(reg(T)) 1 1 M 4 1 1 L „ (n) 1 1 V" , Ai 



x,At I 



Similarly, thanks to further consistency results, the two last terms in (48) can be 
rewritten as 

cT 



-^app 



ip+ / s i </j+r-p(size(T),Ai)C(reg(X))('||/ a pp|| I ,i(Q)+||si|| i i((o i T)xr J v)- 

Gathering the above calculations, we end up with the weak form of the discrete 
equation: 

(49) j^(-^ A *(-)5^ + M,(.)V|, AtU ,r At (.)- V¥> + fc(«*' At (-))v) 

va<p(0,-)+ // I app (p+ / / Si ip 



JJq Jo Jr N 

+ Cr y (size(T), At)(l + || V^ At ul At \\ Ll{Ql 

the constant C depends on the data of the problem, O and reg(T). The second 
equation of the system is analogous, with ' At replaced by uJ ,At and with signs 
changed accordingly. 

Step 5 (Convergences via compactness). All the convergences below are along a 
sub-sequence of a sequence At m , T m of time steps and meshes with size(T m ) + At rn 
tending to zero as m — > oo. In what follows, we drop the subscripts "m" in the 
notation; indeed, after the identification of the limits, the uniqueness of a solution 
to Problem ([I]) , ^ , ^ will permit to suppress the extraction argument. 

From ( |47| and the compactness results in Sections 5.4 5.5 we readily find that 

(50) <; A *(-) -> Ui,e, (V g r xAt <; A *)(-)-^ Vu i>e weakly in L 2 (Q), 

as size(T), At ->■ 0; and u ilB - &, e € L 2 (0, T; V). Because v T < At = uf' At - u^ At (-), 
analogous convergences hold for v' x ' At and its discrete gradient, the corresponding 
limits being v := Ui — u e and Vw, respectively. 



Moreover, if = <9f2, g ije = 0, we can use the compactness result of Section 5.6 



and infer the strong convergence of u T ' At (-) in L^Q); by the preceding remark 
the limit is identified with v. In the case of other boundary conditions, we use the 
local version of Proposition [5~4] (shown in [5] for traditional finite volume schemes; 
the adaptation to DDFV schemes is straightforward) . Up to now we only have the 
i 1 (0, T; Lj oc (Q)) convergence of w x ' At (-) to v. In both cases, the uniform up-to-the- 
boundary estimate ( 46 ) and the interpolation argument yield the strong convergence 
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of v T ' At (-) to v in L r e (Q), for all e > 0, and the weak L r (Q) convergence. In 
particular, thanks to the growth assumption ^ on h we have 

(51) /i(v T ' At (-)) -> h(v), strongly in L\Q) as size(X), At -> 0. 

Step 6 (Passage to the limit in the continuous weak formulation) . In view of the 



properties (501,(51 ) of Step 5, the passage to the limit in (49) and the corresponding 
equation for uJ'^Ms straightforward. We conclude that the limit triple (v,i,u e ,v) 

of (u'* At ,u*' At ,v'*' At ^ is a weak solution of Problem Q,([2]),([3]). In view of the 
uniqueness of a weak solution, we can bypass the "extraction of a sub-sequence" 
part in Step 5. This ends the convergence proof. 

Step 7 (Strong convergences). We will prove that the functions uf ' e ' and their 
discrete gradients converge strongly to u,- ie , Vttj ie , respectively, in L 2 (Q), while 
v *,At converges strongly to v in L r (Q). To this end, we will utilise monotonicity 
arguments to improve the weak convergences to the strong ones. 

By the established weak convergences and the strong L 2 convergence of 
and (for version (C)) of v™ 0,0 to vq, we get 



lim 

3(x),At^0\2 



(52) 



N 



A? 



app > v 



n=0 



, j.// 1,1+1 9x,n+l 



'app 1 



/V 



n=0 



o Jr, 



(sjMi + S e U e ). 



,T,n+l flx.ra+l 



First, as in Step 2, take the discrete solutions U* ^ A< as test functions in the discrete 
equations and subtract the resulting identities. Next, with the help of the regular- 



isation Lemma 2.4 take test functions in the two equations of the system, 



and subtract the resulting identities. Comparing the two relations with the help of 



(52 1, using in addition inequality (41), we infer 



lim I - / \v 

-ize(x),At->0 \ 2 ./si 



< 



M (T) + JJ h(v^ At (-))v*< At (-) 

JJ (m 4 (.)( v-, At ^' At )(-) • ( Vj, At ^' At )(.) 



- m b (.) ( v* >At uf At ) (•) • ( v^ iAt «r At ) (•) 



h(v) v 



Mi VUi ■ VUi + M e Vw e ■ Vlt< 



Furthermore, let us assume for simplicity that L = 0, I = (which means that 
h(0) = and h(r)r > 0, so that the Fatou lemma can be used); to treat the general 
case, use the test function Q(t) := exp(2L(T — t))]lr 0)T )(t) in order to absorb the 

terms containing £|v| 2 C into the term \dt(- 

By the Fatou lemma and properties of weak convergence (with respect to weighted 
vector- valued L 2 (Q) spaces with weights the matrices Mj je > 0), we conclude that 
the above inequality is actually an equality. Using the fact that weak convergence 
plus convergence of norms yields strong convergence in uniformly convex Banach 
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spaces, using an easy refinement of the Fatou lemma]^] (convergence of the integrals 
implies the strong convergence), we conclude that 

(V^, At <; A *)(-) ->■ Vu ije strongly in L 2 (Q), 

h(v TAt (-))v^' At (-) -> h(v)v strongly in L r (Q). 

Using the lower bound in ^ and the Vitali theorem, we infer that ||v T,A *(-)llL r (Q) 
converges to IMIl^q), thus the weak L r (Q) convergence of v' x ' At (-) to v is upgraded 
to the claimed strong convergence. 

Finally, the strong L 2 (Q) convergence of the discrete gradients of u^' At ensures 
a uniform estimate on their translates in time: 

/ f | V T uf; At (t+r, x) - V x M^ At (t, x)\ 2 dxdt -> as r^O, uniformly in % At. 
Jo Jn 

Then the discrete Poincare inequality yields a uniform control of the L 2 (Q) time 
translates 

l-T-T I- 

\ U ie (* + T ' X ) ~ U ie t {t,x)\ 2 dxdt 



o Jn 

of uff * (here we also use the uniform time translates of the discrete Dirichlet data 
<7j X g A *; as usual, the case = is treated separately). Because we can control 
the space translates of uJ' At through the uniform L 2 (Q) estimate of W^uf'^ (see 
5.3 1, we conclude that u*' At converge strongly in L 2 (Q) to to Ui, e . 



Section 



Remark 6.1. Under some stronger proportionality assumptions on the meshes, 
consistency properties similar to those of Proposition |5.1[ i) , (ii) hold not only for test 
functions, but also for functions in L 2 (0, T; iJ 1 (fi)) (cf. [5]). Using the argument 
of [6], we conclude that the discrete solution uf' e converges strongly in L 2 (il) to 
Uie- Indeed, from the discrete Poincare inequality we derive the estimate 



X,At Y7X ml „ cjAt 



\i i > ^ _ pan" Si",,. 
\ U i,e r ° S U he\\ L 2(Q) 



<C(reg(30,n)||V*,A t u^ 1 - V*,A t P* o§^, e || i2(Q) 

and analogous estimates on the meshes SDT* and (for version (C)) 9Ji° . Then the 
consistency and strong convergence of the discrete gradients imply the desired re- 
sult: we find ,A * — u i,e|| L 2(Q) — » as size(T), At — > 0, and so forth. 

6.2. A linearised implicit scheme and its convergence. We follow step by 
step the preceding proof and indicate the modifications needed to take into ac- 
count the linearised-implicit treatment of the ionic current term. We notice that 
throughout the calculations of the preceding proof, 

(53) fr(v T ' At (-)) should be replaced by b{v'^ At {- - At))v T ' At (-); 

where we have set v T,At (t, •) := u T,0 (-) for t £ (— At, 0], and by v x,At (- — At) we 
mean the function (t, x) e Q n- v T ' At (t — At, x). 

Step 1. We cannot get the continuous dependence with the same technique, 
but by induction, we get uniqueness. Indeed, as soon as the uniqueness of v %,n is 
justified, we have 

(6(«*'"(.))t>* ,n+1 (0 -b(v % ' n (-))v % ' n+1 (-))(v % ' n+1 (-) -v T ' n+1 (-)) >0. 

n 



This inequality plays the same role as the non- negativity of the second term in ( 44 1 



2 This result is sometimes referred to as the Schaeffe lemma, and it can be stated as follows: 
fn > 0, f n — > / a.e. on f2, J f„ — > J f as n — > oo j => ^ f n — > / in L 1 (f2) as n — > oo J . 
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Steps 2 and 4. The arguments are unchanged, except for (53) 



Step 3. The estimates (45 1 and (|47[) remain true. Notice that the function b is 




non-negative. Therefore the estimate mm is replaced by the following one: 



(54) 



6(«*. A *(.-At))|t;*" at (0r<C. 



x. At/ 



Now we u se n ew arguments. Namely the discrete Sobolev embedding inequality 
of Section 5.3 yields a uniform L 2 (0, T; L 6 (f2)) bound on v % ' ; then interpolation 
with the L°°((0,T),i 2 (O)) bound (g5j ensures that 

lU.x.At/ 



(55) 



'(■ 



li 10 / 3 (Q) 



< c. 



Step 5. The main difference is the way we ensure the strong L 1 (Q) convergence 
of v T - At and the weak L l (Q) convergence of the associated ionic current term 



(56) 



,x,At / 



(■) :=b(v*> At (.-At))v*^(.) 



%,Ati 



Lv 



T,At 



(■) - 1. 



It is sufficient to treat the nonlinear part of h' T ' (■); thus we can "forget" about 



the two last terms in ( 56 ) . Let us first notice that the definition of b and the growth 
bound ^ on h imply that for some constant j3 — f3(a, L, I) we have 

Hz) < /3(i + izr 2 ). 

Then the assumption r — 2 < 16/3 — 2 = 10/3 made in Theorem |4.5[ii) and the 



uniform L W / 3 (Q) bound (p>5| on v' T ' At ensure the equi-integrability of the functions 



— At)) on Q. Now for any measurable set E C Q, for all S > 0, 



\b(v^(-~At))v^ At (-)\< 



6(w T < At (--Ai)) 



b(v %At {-- At, X ))\V* < A '(-) 



,At I 



Thus estimate (54) and the aforementioned equi-integrability of b(v T ' At (- — At)) 
ensure the equi-integrability of the ionic current term /i' r ' At (-). In particular, from 



(26) we infer L X (Q) bounds on the components of the discrete function h % ' At : 



N 

E 

n=0 



(» 



At \\h 



I 



|L l (n) 



an* ,n+l i 



|L 1 (f2) 



for version (C) , th e term on 9Jf is also controlled. At this stage, the full strength 
of Proposition 5.4 is put into service: indeed, we only have the L 1 control of the 



right-hand side of the discrete evolution equations (22). We infer the strong L X (Q) 
convergence (along a sub-sequence) for v' x ' At . Then from the Vitali theorem and 



the fact that v*> At (- - At) 
convergence of h T ' At (-) to h(v) 



,t,At 



(•) — > in L 1 (Q) and a.e., we get the strong 



Steps 6 and 7. The arguments remain unchanged, taking into account (53). 



7. Numerical experiments 

For the numerical simulations, the bidomain problem ([T]) is reformulated in terms 
of v and u e only; the elimination of ut, thanks to the relation v = Ui~ u e , decreases 
the number of unknowns per primal/dual volume from three to two. 

In terms of v and u e , the parabolic type problem ([T]) is turned into the following 
elliptic-parabolic problem: 




div(M e (x) + M l (x))Vu e +divM,(x)Vu = (t,x) e Q, 
ed t v + e 2 divM e (x)Vu e + h[v] = / ap p (t,x) € Q. 
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For e = 1, problem (57) is equivalent to the original problem ([!]); indeed, the first 
line in ( 57 ) results from the summation of the two lines in ([I ) , with m = v 



The bidomain problem will be considered under formulation ( 57J) throughout this 
section. We point out that numerical schemes associated with formulation (57) 
are equivalent with numerical schemes for the formulation ([I]), following the same 
algebraic operations on the discrete equations. 

A scaling parameter e has also been introduced in ( |57| . Its presence clearly 
makes no difference for the mathematical study of the previous sections, but it 



greatly helps the solutions of (57) to behave as excitation potential waves (which 
waves ( 57 ) is supposed to model) . More precisely, following the analysis in [T7] , such 
a scaling parameter together with a cubic shape for h[v] := v(v — l)(v — a) provide 
a simplified model for spreading of excitation in the myocardium; the parameters 
e and a have been set respectively to 1/50 and 0.2. The way excitation waves are 
generated is detailed in Subsection |7.1| 



The convergence of the DDFV space discretisations for the bidomain problem 



has been justified in Theorem 4.5 for two different discretisations (the fully implicit 
in time and the linearised semi-implicit one) of the ionic current term. Here we 
complement these theoretical studies by the numerical experiments on the third 
(and the most important in practise) case of fully explicit in time discretisation of 
the ionic current term. The implementation of the scheme is detailed in Section [7. 2| 
Although the theoretical study of this scheme is made difficult for technical reasons, 
we do observe convergence numerically. 



The convergence result in Theorem 4.5 involves a comparison between the exact 
solution and the discrete solution in the L 2 (Q) norm, the discrete solution being 
interpreted as the weighted sum of two piecewise constant functions (on the primal 
and on the dual cells). The practical computation of this L 2 distance is difficult. 
Namely, a (coarse) numerical solution on a given (coarse) mesh has to be compared 
with a second (reference) numerical solution computed on a reference mesh aimed 
to reproduce the exact solution, unknown in practise. The reference mesh will 
not be here a refinement of the coarse one, thus the precise computation of an 
L 2 norm between the coarse and the reference numerical solutions, thought as 
piecewise constant functions, is an awful task. Therefore we use a slightly different 



convergence indicator, as presented in Section 7.3 The quantity (63) provides 
a more convenient measure of the square of the L 2 (Q) norm. The convergence 
with respect to this variant of discrete L 2 norm will be studied here. Notice that 
under reasonable regularity assumptions on the mesh, the two discrete L 2 norms 
are equivalent. 

Convergence will be studied both in 2D and in 3D. The introduction of the 
two dimensional case is mostly intended to confirm the results in dimension three: 
because of the numerical facilities in dimension two (smaller growth of the problem 
size under refinement), it allows a deeper insight into the asymptotic behaviours. 



7.1. Settings. Problem (57) is considered for a reaction term v i-> h[v] set on the 
domain SI = [0, l] d , d = 2,3 denoting the space dimension. We consider solutions 
under the form of excitation potential waves spreading from the domain centre 
towards its boundary, as depicted in Figure|3]in the two dimensional case, relatively 
to some medium anisotropy detailed below. 

Excitation is initiated by applying a centred stimulation during a short period of 
time, precisely: I app (x, t) — 0.9 for 1 < t < 1.1 and \x — xq\ < 0.1 (x denoting the 
centre of Q) and I app (x, t) = otherwise. The initial condition for v is uniformly set 
to 0. A homogeneous Neumann boundary condition is considered on dtt, uniqueness 
is ensured by adding the normalisation condition Q on u e . The domain Q is 
assumed to be composed of a bundle of parallel horizontal muscular fibres, resulting 
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Figure 3. Reference solution in the two-dimensional case. 
Above: spreading of the transmembrane potential v excitation 
wave. From left to right, the three pictures correspond to times 
i=0.2, 0.6 and 1.2 after the stimulation initiation, stimulation du- 
ration being 0.1. Below: associated extracellular potential u e . 



in the following choice for the anisotropy tensors Mj(i) and M e (x): 

(58) Mi(a:) = Diag(A<, A*, A*) , M,(i) = Diag(Ai, A*, A*) , 

the values for the longitudinal (I) and transverse (t) conductivities for the intra and 
extra-cellular medias have been taken from |42) : the resulting anisotropy ratios 
for the intra and extra-cellular medias respectively are 9.0 and 2.0 between the 
longitudinal and transverse directions. 

The numerical solution for the transmembrane potential v takes the form of an 
excitation wave propagating across the domain from the stimulation site towards 
the boundary and from the rest potential v — to the activation potential v = 1. A 
sharp but smooth wavefront for v displays an elliptic shape away from the boundary, 
which is induced by the media anisotropy. A reference solution is generated on a 
mesh using a 1 147 933 (resp. 479 873) nodes in dimension 3 (resp. 2). 

7.2. Implementation. Let us fix a mesh %. For simplicity, discrete functions 
w,v £ R T will also be considered as one-column real matrices in this subsection, 
w T , v T denoting their transpose one-row real matrices. Let us first introduce the 
mass matrix (diagonal here) A £ Mat(M T ): 



\fw, v £ 



w, v 



= w Av. 



Relative to (58 1, uniform discrete tensors Mf e are considered here, with value 
Mf = Diag(A<, A*, A*) , M* = Diag(Ai, A*, A*) , 

on each diamond d. The discrete gradient being defined relative to the homoge- 
neous Neumann boundary condition, and simply denoted by V x , the two stiffness 
matrices S, and S P are introduced as: 



These stiffness matrices are positive, symmetric matrices, although not definite 
since a Neumann homogeneous boundary condition is considered. 
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The following semi-implicit Euler scheme is considered: given u n ,J" pp £ 



determine v n+1 ,u" +1 £ K -1 such that, 



(59) 



div x [(M e r +Mf)V*«* +1 ] +div T [Mf V T w" +1 ] =0, 



,,n+l 



At 



e 2 div X [MJ V T < +1 ] + h[v n ] = I, 



app- 



12 



of each line in ( 59 1 with all test func- 



Writing separately the scalar product 

tions w £ K -1 and using the discrete duality ( |11[ ) leads to the following equivalent 
formulation written in matrix form: 



,71+1 







— £AiS e 



A 



(60) M = , M := 

A[v n + At(i: pp - h[v n ])/e] 

To ensure uniqueness on u e , the normalisation condition Q is discretised as: 



(61) 



,ri+l 



= = 



n+l 



where U™° (resp. U m * ) is the discrete function equal to 1 relatively to each primal 
(resp. dual) control volumes and to elsewhere. 

The implementation of ( 60 1 therefore reads as a three-step algorithm: 

at each time step, 

(1) compute y2 := A[v n + Ai(I" — h[v n ])/e]. The matrix A being diagonal, 
computations for this step are cheap; 

(2) determine a solution x — (u" +1 ,v n+1 ) T to the global system Mx = y for 

y-(o,2/ 2 ) T 



(3) normalise using condition (61). For the same reason as for step one, 

this step is a cheap one. 

Because of the large size of the considered problem (1.1 million of nodes in 3D 
for the most refined mesh, i.e. 2.2 million of lines for the matrix M) and because 
of the relatively non-compact sparsity pattern for M in 3D, step 2 is not an easy 
task. Therefore, a careful attention has to be paid to the preconditioning of M: 
the strategy adopted here is detailed in [50] . 

7.3. Numerical tests and results. The convergence of the DDFV scheme is nu- 



merically analysed comparing the reference solution described in Subsection 7.1 
with numerical solutions obtained on coarser meshes. In 3D, four tetrahedral 
meshes have been considered: from 2 559 to 1 147 933 nodes, between two meshes 
the mesh size is divided by 2, two successive meshes are not obtained via refine- 
ment. In 2D, six meshes are used: from 489 to 479 873 nodes. The time step At 
is also divided by two each time the space resolution is divided by 2; the starting 
time step (on the coarsest mesh) is 0.02. 

To compare numerical solutions defined on different meshes, a projection is 
needed: this is done as follows. Let T r and T c be the reference mesh and a coarser 
mesh respectively, and let R Tr , M Xc respectively denote the associated spaces of 
discrete functions. Consider the simplicial mesh S r (respectively S c ) whose cells 
are obtained by cutting all diamonds of T r (resp. T c ) in two along the interface. 
A discrete function u r £ M. Tr (resp. u c £ R Tc ) consists in one scalar associated to 
each vertex and each cell centre of the mesh X 1 " (resp. T c ), thus to each vertex of 
S r (resp. S c ). It is therefore natural to associate to u c (resp. u r ) the continuous 
function u r (resp. u c ) piecewise affine on the cells of S r (resp. S c ) and whose 
values at the vertices of S r (resp. S c ) are given by the discrete function u r (resp. 
u c ). A projection u c ^ r £ K Tr of a (coarse) discrete function u c £ M Xc is then 
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# nodes 


errors en, 2 


489 


1.611 


1913 


9.130 10~ 2 


7 569 


1.776 10~ 2 


30 113 


4.850 10- 3 


120 129 


1.139 10- 3 


479 873 


reference 



# nodes 


errors e , 2 


2 559 


1.110 


19 500 


8.195 l(T a 


148 242 


1.281 10~ 2 


1 147 933 


reference 



2D case 3D case 

Table 1. Activation time mappings convergence. The errors are 



relative errors in L (ft) norm as denned in (62) 



simply defined by computing the values of the function u c on the vertices of S r . 
The relative error between u r € M Tr and u c 6 R Xc in L 2 (il) norm is defined as 



■ ^ _ fa 



Numerically, these integrals are evaluated using an order two Gauss quadrature on 
the cells of S r , leading to an exact evaluation up to rounding errors. 
The three following tests have been performed. 




Figure 4. Activation time in dimension 2 for three different 
meshes, the isolines (in black) are separated by 1/3 unit of time. 
The stimulation is initiated at time t = 1. From the left (coarsest 
mesh) to the right (reference solution) the three different activa- 
tion time mappings have been computed on meshes with 439, 7569 
and 479 873 nodes respectively. 

7.3.1. Test 1 ; activation time convergence. The activation time mapping <j> : Sl^I 
is defined at each point x as the time cj>(x) = t such that the transmembrane 
potential v(x,t) = s for the threshold value s := 0.9. The value 4>(x) tells us at 
what time the excitation wave reaches the point x, the activation time mapping 
thus is of crucial importance in terms of physiological interpretation of the model. 
Activation time in 2D computed on various meshes are depicted on Figure |3J The 
discrepancy between the activation mappings (jf and <\f computed at the reference 
and coarse levels respectively is evaluated using the relative error in the I?(SX) 



norm defined in ( 62 ) . Numerical results for activation time convergence are given 
in Table [3 

Convergence is numerically observed here both in 2D and in 3D. Moreover, the 
figures obtained in the two dimensional case indicate an order two convergence 
relatively to the mesh size. Such a conclusion, although plausible, is not possible in 
the three dimensional case: to be observed it would require a much finer reference 
mesh which is not affordable in terms of computational effort. 
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# nodes 


f=1.2 £ = 1.6 £ = 2.2 


f = 1.2 £ = 1.6 £ = 2.2 


489 


0.51 0.37 0.45 


0.45 0.34 0.59 


1913 


0.24 8.87 10~ 2 0.14 


0.22 8.27 10~ 2 0.17 


7569 


0.13 6.01 10" 2 1.42 10" 2 


0.12 5.38 10' 2 1.99 10" 2 


30113 


6.24 10~ 2 3.29 10~ 2 1.28 10~ 2 


5.62 10" 2 2.89 10" 2 1.71 10~ 2 


120129 


1.25 10~ 2 5.26 10~ 3 1.74 10" 3 


1.13 10~ 2 5.31 10" 3 2.86 10~ 3 



2D case 



# nodes 


£ = 1.2 £=1.6 £ = 2.2 


£=1.2 £ = 1.6 £ = 2.2 


2559 


0.43 0.57 0.65 


0.42 0.63 0.91 


19 500 


0.22 0.14 0.20 


0.19 0.17 0.28 


148 242 


0.10 6.21 10~ 2 2.88 10" 2 


9.45 10" 2 6.40 10" 2 4.33 10~ 2 



3D case 

Table 2. Convergence of the transmembrane potential v and of 
the extracellular potential u e at three fixed times: £ = 1.2, 1.6 and 
2.2. The reference solution for these chosen times are depicted in 
Figure [3] Above (resp. below) are reported the errors in the two 
(resp. three) dimensional case. On each table line the three first 
figures after the number of nodes correspond to the errors on v, 
whereas the three last ones correspond to u e . Errors are relative 



errors in L (f2) norm as defined in (63) 



7.3.2. Test 2; space convergence. Let us denote by v r and u r e , (resp. v c and u^) 
the transmembrane potential and extracellular potential computed on the reference 
mesh (resp. a coarse mesh). The discrepancy between v c , v r and u r e , u% at a chosen 
time £ has been computed using the relative error in the L 2 (fl) norm (62). Three 
fixed times £ have been considered: t = 1.2, 1.6 and 2.2, corresponding to the 
reference solution depicted in Figure [3] 

Because of the particular wavefront-like shape of the solution, this error can 
be geometrically reinterpreted as follows. Consider at time £ the sub-region of ft 
that is activated according to the reference solution u r e but not activated according 
to the coarse solution u%. The numerator in (63) simply measures the square 
root of the area of this sub-region. Using the elliptic shape of activated regions, 
one gets that eQ^{u r ,u c ) measures the square root of the relative error on the 
wavefront propagation velocity (more precisely the square root of the sum of the 
axial and transverse wavefront propagation velocities relative errors). This error, as 
in test case 1, is of prime physiological importance. Numerical results for this test 
are displayed in Table [2] Although convergence is well illustrated, no particular 
asymptotic behaviour can be inferred from these results. 



7.3.3. Test 3; space and time convergence. A numerical space and time convergence 
indicator eQ.2 is introduced here, aiming to reproduce an L 2 (Q) relative error be- 
tween a coarse and the reference solution (Q = (0, T) x fi). Convergence is mea- 
sured using this indicator, and this third test therefore is intended to numerically 
illustrate the convergence result of Theorem |4.5| 

The transmembrane potentials v r and v c have been recorded at the same times, 
namely £„ = n8t, for St = 1/100 unit of time, n = 0, . . . ,T/6t and T = 3. The 
corresponding numerical solutions are denoted v r n and v°. A relative error in the 
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# nodes 


errors e Q:2 


489 


0.481 


1913 


0.237 


7569 


6.469 10" 2 


30113 


1.746 10~ 2 


120129 


4.167 10- 3 


479 873 


reference 



# nodes 


errors e Q:2 


2559 


0.673 


19 500 


0.219 


148 242 


4.920 10~ 2 


1 147 933 


reference 



2D case 3D case 

Table 3. Space and time convergence for the transmembrane po- 
tential v. Errors are relative errors in the L 2 (Q) norm as defined 
in (|63l 



L 2 (Q) norm between v r and v c is introduced as follows: 

/ co \ / r c\2 En=0 In \ V n — V n \ 2 dxSt 

(63) e Qt2 (v,v):= „jv } , - r , 2 , x . - N = T/St. 

Numerical results are given in Table [3} convergence both in 2D and in 3D is 
observed. The two dimensional case results indicate an order two convergence 
relatively to the mesh size. As in the first test case, such a conclusion is not 
possible in the three dimensional case: a much finer reference mesh (not affordable) 
would be needed. 
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